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0^ ' We discuss various aspects of the post-Newtonian approximation in general relativity. 

After presenting the foundation based on the Newtonian limit, we use the (3+1) formalism to 
formulate the post-Newtonian approximation for the perfect fluid. As an application we show 
the method for constructing the equilibrium configuration of nonaxisymmetric uniformly 
rotating fluid. We also discuss the gravitational waves including tail from post-Newtonian 
systems. 



§1. Introduction 



00 ' The motion and associated emission of gravitational waves (GW) of self gravi- 

■ fating systems liave been a main researcli interest in general relativity. The problem 

\^ ', is complicated conceptually as well as mathematically because of the nonlinearity 

' of Einstein's equations. There is no hope in any foreseeable future to have exact 



solutions describing motions of arbitrary shaped, massive bodies so that we have to 
adopt some sort of approximation scheme for solving Einstein's equations to study 
such problems. In the past years many types of approximation schemes have been 
^ developed depending on the nature of the system under consideration. Here we 

t^', shall focus on a particular scheme called the post-Newtonian (PN) approximation. 

There are many systems in astrophysics where Newtonian gravity is dominant, but 
general relativistic gravity plays also important roles in their evolution. For such 
systems it would be nice to have an approximation scheme which gives Newtonian 
' description in the lowest order and general relativistic effects as higher order per- 

turbations. The post-Newtonian approximation is perfectly suited for this purpose. 
HistoricallvEinstein computed first the post-Newtonian effects, the precession of the 
perihelion Lj), but a systematic study of the post-Newtonian approximaiioji^was iiDt 
made until the series of papers by Chandrasekhar and associatescil'E3)'c3)'EfP'E3. 
Now it is widely known that the post-Newtonian approximation is important in an- 
alyzinga number of relativistic problems, such as the equaiipnpf motion of binary 
pulsar Ej)'c2P' 13 'Ej) solar-system tests of general relativityl^'tS*, and gravitational 
radiation reaction@'0. 

Any approximation scheme nesseciates a small parameter(s) characterizing the 
nature of the system under consideration. Typical parameter which most of schemes 
adopt is the magnitude of the metric deviation away from a certain background 
metric. In particular if the background is Minkowski spacetime and there is no other 
parameter, the scheme is sometimes called as the post-Minkowskian approximation 
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in the sense that the constructed spacetime reduces to Minkowski spacetime in the 
limit as the parameter tends to zero. This hmit is called as the weak field limit. 
In the case of the post-Newtonian approximation the background spacetime is also 
Minkowski spacetime, but there is another small parameter, that is, the typical 
velocity of the system divided by the speed of light which we call e henceforth. 
These two parameters (the deviation away from the flat metric and the velocity) 
have to have a certain relation in the following sense. As the gravitational field 
gets weaker, all velocities and forces characteristic of the material systems become 
smaller, in order to permit the weakening gravity to remain an important effect on 
the system's dynamics. For example in case of a binary system, the typical velocity 
would be the orbital velocity v/c ^ e and the deviation from the flat metric would 
be the Newtonian potential, say <P. Then these are related by <P/c^ ~ jt? ~ 
which guarantees that the system is bounded by its own gravity. 

In the post-Newtonian approximation, the equations in general relativity take 
the form of Newton's equations in an appropriate limit as e ^ 0. Such a limit is called 
as the Newtonian limit and it will be the bases of constructing the post-Newtonian 
approximation. However the limit is not in any sense trivial since the limit may be 
thought of two limits tied together as just described. It is also worth noting that 
the Newtonian limit cannot be uniform everywhere for all time. For example any 
compact binary systems, no matter how weak the gravity between components and 
slow its orbital motion, will eventually spiral together due to backreaction of the 
emission of gravitational waves. As the result the effects of its Newtonian gravity 
will be swamped by those of its gravitational waves. This will mean that higher 
order effect of the post-Newtonian approximation eventually dominates the lowest 
order Newtonian dynamics and thus if the post-Newtonian approximation is not 
carefully xipnstructed, this effect can lead to many formal problems, such as divergent 
integrals It has been shown that such divergences may be avoided by carefully 
defining Newtonian limited'. Moreover, the use of such limit provides us a strong 
indication that the post-Newtonian hierarchy is an asymptotic approximation to 
general relativity cIP . Therefore we shall first discuss in this paper the Newtonian 
limit and how to construct the post-Newtonian hierarchy before attacking practical 
problems in later sections. 

Before going into the details, we mention the reason for the growth of interest 
in the post-Newtonian approximation in recent years. Certainly the discovery of 
the binary neutron star system PSR 1913+16 was a strong reason to have renewed 
interest in the post-Newtonian approximation since it is the first svstem in which 
general relativitistic gravity plays fundamental roles in its evolution Particularly 
indirect discovery of gravitational wave by the observation of the period shortening 



led to many fruitful 
emission in 1980'sizl)'l3) 



Rs of the equation of motion with gravitational radiation 
y . The effect of radiation reaction appears as the form of 
the potential force at the order of e' higher than Newtonian force in the equation 
of motion. There have been various attempts to^how the va^^ ilf tJie ^-called 
quadrupole formula for the radiation reaction B)'C§'l^'E^'E3)'£ZP'0'Lil'O*'^ 

At that time, however, no serious attempts had been made for the study of 
higher order effects in the equation of motion. The situation changed gradually in 
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1990's because of the increasing expectation of direct detection of gravitational waves 
by Kiloin£ter-size interferometric gravitational wave detectors, such as LIGOUP, 
VIRGOII3) and TAMA tip now under construction. Coalescing binary neutron stars 
are the most promising candidate of sources of gravitational waves for such detec- 
tors. The reasons are that (1) we expect tcLdetect the signal of coalescence of binary 
neutron stars about several times per yearll^), and (2) the waveform from coalescing 
binaries can be predicted with a high accuracy compared with other sources 
Informations carried Jay^'avitational waves tell us not only various nhysical param- 
eters of neutron starty)'ElP, but also the cosmological parameters e§)'H'eHP'E3) if and 
only if we can make a detailed comparison between the observed signal with theoret- 
ical prediction during the epoch of the so-called inspiraling phase where the orbital 
separation is much larger than the radius of component stars This is the place 
where the post-Newtonian approximation may be applied to make a theoretical tem- 
plate for gravitational waves. The problem is that in order to make any meaningful 
comparison between theory and observation we need to know the detailed waveforms 
generated by the motion up to 4PN order which is of order higher than the Newto- 
nian order E3^. This request from gravitational wave astronomy forces us to construct 
higher order post-Newtonian approximation. 

It should be mentioned that Blanchet and Damour have developed a system- 
atic scheme to calculate the waveform at higher order where the post-Minkowskian 
approximation is used to construct the external field and the post-jMewtx^ i aij 



proximation is used to construct the field near the material source^' 
Blanchet has obtained the waveform unto the 2.5 PN order which is of order 
higher than the lowest qii dra upole wavellit'i'*'& by using the equation of motion up 
to that orderE2)'E2l''E9)'L2P'E2P. We shall also study the waveform from different point 
of view because of its importance. We obtain essentially the same results including 
the tail term with them. 

There is another phase in the evolution of the binary neutron star which we can 
withdraw useful informations from the observed signals of gravitational waves. That 
is the intermediate phase between the inspiralling phase and the event of merging. 
In the inspiralling phase the components are usually treated as point particles, while 
the extendedness becomes important in the intermediate phase and the merging. 
Pull general relativistic hydrodynamic simulation is needed for the understanding of 
the merging. The material initial data for the simulation will be the density and the 
velocity distribution of the fluid. The post-Newtonian approximation will be used to 
construct the data. Furthermore the post-Newtonian approximation is able to take 
into account of the tidal effect which is important for the orbital instability. Recently, 
Lai, Rasio and Shapiro l3^'E5) have pointed out that such a tidal coupling of binary 
neutron stars is very important for their evolution in the final merging phase because 
the tidal effect causes the instability of the circular motion even in the Newtonian 
gravity. Also important is the general relativistic gravity because in the intermediate 
phase, the orbital separation is about ten times as small as the Schwarzschild radius 
of the system. Thus, we need not only a hydrodynamic treatment, but also general 
relativistic one in order to study the final stage of the evolution of binary neutron 
stars. Since the usual formulation of the post-Newtonian approximation is based on 
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the particular gauge condition which is not so convenient for numerical purposes, 
we shall reformulate it by the (3+1) formalism often used in numerical approach of 
general relativity. 

This paper is organized as follows. In section 2 we introduce the Newtonian 
limit in general relativity and present how to construct the post-Newtonian hierar- 
chy. There we mention how to avoid divergent integrals which appears at higher 
order in the previous treatments, and we also discuss how to incorporate strong in- 
ternal gravity in the post-Newtonian approximation. Next we reformulate the post- 
Newtonian approximation appropriate for numerical treatment in section 3. There 
we shall adopt the (3+1) formalism which is frequently used in numerical relativity. 
Based on the formalism developed in section 3, we present a formulation for con- 
structing numerically equilibrium solutions of uniformly rotating fluid in the 2PN 
approximation in section 4. We shall also discuss the propagation of gravitational 
waves from slow motion systems in section 5. 



§2. Foundation of the post-Newtonian approximation 



Since the Newtonian limit is the basis of the post-Newtonian approximation, we 
shall first formulate the Newtonian limit. We shall follow the formulation by Futa- 
mase and SchutzElP. We will not mention other formulation of Newtonian limit by 
Ehlers§)S, because it has not yet used to construct the post-Newtonian approxi- 
mation. 

2.1. Newtonian limit along a regular asymptotic Newtonian sequence 

The formulation is based on the observation that any asymptotic approximation 
of any theory need a sequence of solutions of the basic equations of the theorycZP'EHP . 
Namely, if we write the equations in abstract form as 

E{g) = 0, (2-1) 

for an unknown function g, one would like to have a one-parameter (or possibly 
multiparameter) family of solutions, 

E{g{X)) = 0, (2-2) 

where A is some parameter. Asymptotic approximation then says that function /(A) 
approximates g{X) to order A^ if |/(A) — g{X)\/X'P ^ as A — > 0. We choose the 
sequence of solutions with appropriate properties in such a way that the properties 
reflect the character of the system under consideration. 

We shall formulate the post-Newtonian approximation according to the general 
idea just described. As stated in the introduction, we would like to have an approx- 
imation which applies the systems whose motion is described almost by Newtonian 
theory. Thus we need a sequence of solutions of Einstein's equations parameterized 
by e (the typical velocity of the system divided by the speed of light) which has the 
Newtonian character as e — > 0. 

The Newtonian character is most conveniently described by the following scal- 
ing law. The Newtonian equations involve six variables, density (p), pressure (P), 
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gravitational potential (^), and velocity {v^,i 



1,2,3): 



- Anp = 0, 



(2-3) 



dtp + Vi{pv') = 0, 
pdtv' + pv^Vjv' + V'P + pV'^ = 0, 



(2-4) 
(2-5) 



supplemented by an equation of state. For simplicity we have considered perfect 
fluid. Here we have set G = 1. 

It can be seen that the variables t), t), t), v'^{x^ ,t)} obeying the 

above equations satisfy the following scaling law. 



One can easily understand the meaning of this scaling by noticing that e is the 
magnitude of typical velocity (divided by the speed of light). Then the magnitude of 
the gravitational potential will be of order because of the balance between gravity 
and centrifugal force. The scaling of the time variable expresses the fact that the 
weaker the gravity (e — 0) the longer the time scale. 

Thus we wish to have a sequence of solutions of Einstein's equations which has 
the above scaling as e ^ 0. We shall also take the point of view that the sequence 
of solutions is determined by the appropriate sequence of initial data. This has a 
practical advantage because there will be no solution of Einstein's equations which 
satisfy the above scaling (2-6) even as e ^ 0. This is because Einstein's equations 
are nonlinear in the field variables so it will not be possible to enforce the scaling 
everywhere in spacetime. We shall therefore impose it only on the initial data for 
the solution of the sequence. 

Here we shall first make a general discussion on the formulation of the post- 
Newtonian approximation independent of any initial value formalism and then present 
the concrete treatment in the harmonic coordinate. The condition is used because of 
its popularity and some advantages in the generalization to the systems with strong 
internal gravity. 

As the initial data for the matter we shall take the same data set in Newtonian 
case, namely, the density p, the pressure P, and the coordinate velocity v^. In most 
of the application, we usually assume a simple equation of state which relates the 
pressure to the density. The initial data for the gravitational field are Qf^uidgf^u/dt. 
Since general relativity is overdetermined system, there will be constraint equations 
among the initial data for the field. We shall write the free data for the field as 
{Qij,Pij) whose explicit forms depend on the coordinate condition one assumes. In 
any coordinate we shall assume these data for the field vanish since we are interested 
in the evolution of an isolated system by its own gravitational interaction. It is 
expected that these choice corresponds to the absence of radiation far away from the 



p{x\t) e^p{x\et), 

P{x\t) ^ €^P{x\et), 
v\x^,t)^ev\x^,et), 
^{x\t) e^^{x\€t). 



(2-6) 
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source. Thus we choose the foUowing initial data which is indicated by Newtonian 
scahng: 

p{t = 0,x\e) = e'^a{x'), 
P{t = 0,x\e) = e%{x'), 
v\t = 0,x'',e) = ec'{x''), 
Q,j{t = 0,x\e) = 0, 

Pij{t = 0,x\e) = 0, (2-7) 

where the functions a,b , and c* are C°° functions that have compact support con- 
tained entirely within a sphere of a finite radius. 

Corresponding to the above data, we have a one-parameter set of spacetime 
parameterized by e. It may be helpful to visualize the set as a fiber bundle, with 
base space R the real line (coordinate e) and fiber the spacetime (coordinates 
t, X*). The fiber e = is Minkowski spacetime since it is defined by zero data. In 
the following we shall assume that the solutions are sufficiently smooth functions 
of e for small e 7^ 0. We wish to take the limit e — > along the sequence. The 
limit is, however, not unique and is defined by giving a smooth nowhere vanishing 
vector field on the fiber bundle which is nowhere tangent to each fiber Elt'EB'. The 
integral curves of the vector field give a correspondence between points in different 
fibers, namely events in different spacetime with different values of e. Remembering 
the Newtonian scaling of the time variable in the limit, we introduce the Newtonian 
dynamical time: 

r = et, (2-8) 

and define the integral curve as the curve on which r and stay constant. In fact 
if we take the limit e — > along this curve, the orbital period of the binary system 
with e = 0.01 is 10 times of that of the system with e = 0.1 as expected from the 
Newtonian scaling. This is what we define as the Newtonian limit. Notice that this 
map never reaches the fiber e = (Minkowski spacetime). There is no pure vacuum 
Newtonian limit as expected. 

In the following we assume the existence of such a sequence of solutions con- 
structed by the initial data satisfying the above scaling with respect to e. We shall 
call such a sequence is a regular asymptotically Newtonian sequence. We have to 
make further mathematical assumptions about the sequence to make explicit calcu- 
lations. We will not go into details of them partly because to prove the assumptions 
we need a deep understanding of the existence and uniqueness properties of Cauchy 
problem of Einstein's equations with perfect fiuids of compact support which are not 
available at present. 

2.2. Post- Newtonian Hierarchy 

We shall now define the Newtonian, post-Newtonian and higher approximations 
of various quantities as the appropriate higher tangents of the corresponding quanti- 
ties to the above integral curve at e = 0. For example the hierarchy of approximations 
for the spacetime metrics can be expressed as follows: 

g^^{e,T,x') = g^^{0,T,x') + e{Cvg,,u){0,T,x') 
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+^e\Clg^,){0,T,x') + + ^^{C'^g^,,){0,T,x') + Rn+i, (2-9) 



where Cy is the Lie derivative with respect to the tangent vector of the curves defined 
above, and the remainder term R^'^i is 



Taylor's theorem guarantees that the series is asymptotic expansion about e = 
under certain assumptions mentioned above. It might be useful to point out that 
the above definition of the approximation scheme may be formulated purely geomet- 
rically in terms of jet bundle. 

The above definition of the post-Newtonian hierarchy gives us asymptotic series 
in which each term in the series is manifestly finite. This is based on the e depen- 
dence of the domain of dependence of the field point (r, x^). The region is finite with 
finite values of e, and the diameter of the region increases as as e — ^ 0. With- 
out this linkage of the region with the expansion parameter e, the post-Newtonian 
approximation leads to divergences in the higher orders. This is closely related with 
the retarded expansion. Namely, it is assumed that the slow motion assumption 
enabled one to Taylor expand the retarded integrals in retarded time such as 



and assign the second term to a higher order because of its explicit e in front. This 
is incorrect because r e^^ as e ^ and thus er is not uniformly small in the 
Newtonian limit. Only if the integrand falls off sufficiently fast, can retardation be 
ignored. This happens in the lower-order PN terms. But at some higher order there 
appear many terms which do not fall off sufficiently fast because of the nonlinearity of 
Einstein's equations. This is the reason that the formal PN approximation produces 
the divergent integrals. It turns out that such divergence appears at 3PN order 
indicating the breakdown of the PN approximation in the harmonic coordinate. 
This sort of divergence may be eliminated if we remember that the upper bound 
of integration does depend on e as e^^. Thus we would get something like e"lne 
instead of e" In oo in the usual approach. This shows that the asymptotic Newtonian 
sequence is not differentiable in e at e = 0, but there are no divergence in the 
expansion and it has still an asymptotic approximation in e that involves logarithms. 

2.3. Explicit calculation in Harmonic coordinate 

Here we shall use the above formalism to make explicit calculation in the har- 
monic coordinate. The reduced Einstein's equation in the harmonic condition is 
written as 




(2-10) 




(2-11) 




(2-12) 
(2-13) 



where 



(2-14) 
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e"(' = {-g){T»P + tl1), (2-15) 

where t^'^ is the Landau-Lifehitz pseudotensor@). In the following we shall choose 
an isentropic perfect fluid for T"'^ which is enough for most of applications. 

T'^P = {p + pn + P)u''u^ + Pg''^, (2-16) 

where p is the rest mass density, 77 the internal energy, P the pressure, the four 
velocity of the fluid with normalization; 

g^pu^u^ = -1. (2-17) 

The conservation of the energy and momentum is expressed as 

V^T"^ = 0. (2-18) 

Defining the gravitational field variable as 

l^u ^ ^^,u _ ^_g^i/2gf.u^ (2-19) 

where r]'^'^ is the Minkowski metric, the reduced Einstein's equation ( |2-12| ) and the 
gauge condition ( 2-13| ) take the following form; 

_ h'"^)h'"'^^^ = -IQ-ne^"" + (2-20) 

h^'^ = 0. (2-21) 

Thus the characteristics is determined by the operator (r/"'^ — h°'^)dadfs, and thus 
the light cone deviates from that in the flat spacetime. We shall use this form of 
the reduced Einstein's equations in the calculation of wave form far away from the 
source because the deviation plays a fundamental role there. However in the study 
of the gravitational field near the source it is not neccesary to consider the deviation 
of the light cone away from the flat one and thus it is convenient to use the following 
form of the reduced Einstein's equationscP. 

r/'^'/i^V = -16^^"^ (2-22) 

where 

^ 0a/3 ^ ^"/^A'^'^^, (2-23) 

^a^3^^u ^ (i6vr)-i(/i°^/i'3'^ - h''^h^"'). (2-24) 
Equations ( p-21| ) and (2-22) together imply the conservation law 

= 0. (2-25) 
We shall take as our variables the set {p, P,v^ ,h°'^}, with the definition 

v' = u'/u^. (2-26) 
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The time component of 4- velocity is determined from Eq.( |2T7] ). To make a well- 
defined system of equations we must add the conservation law for number density n 
which is some function of the density p and pressure P: 

Va(nu") = 0. (2-27) 
Equations ( p-25 ) and ( 2-27 ) imply that the flow is adiabatic. The role of the equation 



of state is played by the arbitrary function n{p,p). 

Initial data for the above set of equations are h'^^ ,h°'^Q, p, P, and v^, but not 
all these data are independent because of the existence of constraint equations. 
Equations. ( |2 -211 ) and ( |2-22 ) imply the four constraint equations among the initial 



Z\/i"0 + 167ryl"°-5'J/i^.'° = 0, (2-28) 



data for the field. 

where A is the Laplacian in the flat space. We shall choose h'^^ and n-'Q as free 
data and solve Eq.(|]2|) for /i"°(a = 0,..,3) and Eq. (|2^) for ,o- Of course 



these constraints cannot be solved explicitly, since contains but they can 
be solved iteratively as explained below. As discussed above, we shall assume that 
the free data h^^ and /I'^'g for the field vanish. One can show that such initial data 
satisfy the O'Murchadha and York criterion for the absence of radiation far away 
from the sourceB^. 

In the actual calculation, it is convenient to use an expression with explicit 
dependence of e. The harmonic condition allows us to have such expression in terms 
of the retarded integral. 

/i^^(e,r,x^) = 4 / d^yA'"'{T-er,y\e)/r + h'ii'{e,T,x'), (2-29) 

where r = — and C(e,r, x*) is the past flat light cone of event (t, x*) in the 
spacetime given by e, truncated where it intersects with the initial hypersurface 
r = 0. h'^ is the unique solution of the homogeneous equation 

= 0. (2-30) 

We shall henceforce ignore the homogeneous solution because they play no important 
roles. Because of the e dependence of the integral region, the domains of integration 
are finite as long as e 7^ and their diameter increases as as e — > 0. 

Given the formal expression (2-29) in terms of initial data (2-7), we can take 



the Lie derivative and evaluate these derivatives at e = 0. The Lie derivative is 
nothing but partial derivative with respect to e in the coordinate system for the 
fiber bundle given by (e,r, x*). Accordingly one should convert all the time indices 
to r indices. For example, T"^"^ = e^T** which is of order since T** p is of order 
e"^. Similarly the other components of stress energy tensor T"^* = eT** and T^^ are of 
order e'^ as well. Thus we expect the first non-vanishing derivative in ( |2-29| ) will be 
the forth-derivatives. In fact we find 

(4)^-(r,x^)=4/ (H)^d3^, (2.31) 
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^ ' Jr3 r 



where we have adopted the notation 



(2-33) 



n)f{r, x^) = A lim r, x*), (2-34) 



and 



(A)ttL = ^ii4)h^^'\4)h'^'' - ls^\4)h^^'\4)h^:k)- (2-35) 



In the above calculation we have taken the point of view that h^'^ is a tensor 
field, defined by giving its components in the assumed harmonic coordinate as the 
difference between the tensor density \/—gg^'^ and rj^'^ . 

The conservation law (2-25) also has its first nonvanishing derivative at this 



order, which are 

{2)p{r, x^)^r + ((2)P(r, x'')^i)v\t, x^))^i = 0, (2-36) 

W(l)V'),r + {(2)P(l)V\l)V')j + i4)P'' - \i2)Pi4)h^^'' = 0- (2-37) 

Equations (2-31), ( |2-36| ), and ( 2-37D consist of Newtonian theory of gravity. Thus the 



lowest non-vanishing derivative with respect to e is indeed Newtonian theory, and 
the IPN and 2PN equations emerge from sixth and eighth derivatives, respectively, 
in the conservation law ( p-25| ). At the next derivative, the quadrupole radiation 
reaction term comes out. 

2.4. Strong point particle limit 

If we wish to apply the post-Newtonian approximation to the inspiraling phase 
of binary neutron stars, the strong internal gravity must be taken into account. The 
usual post-Newtonian approximation assumes explicitly the weakness of gravitational 
field everywhere including inside the material source. It is argued by appealing 
the strong equivalence principle that the external gravitational field which governs 
the orbital motion of the binary system is independent of internal structure of the 
components up to tidal interaction. Thus it is expected that the results obtained 
under the assumption of weak gravity also apply for the case of neutron star binary. 
Experimental evidence for_the strong equivalence principle is obtained only for the 
system with weak gravity B3)'L3), but no experiment is available in case with strong 
internal gravity at present. 

In theoretical aspect, the theory of extended object in general relativity!^ is still 
in preliminary stage for the application to realistic systems. Matched asymptotic 
expansion techninue has been used to treat the system with strong gravity in certain 
situations Another way to handle with strong internal gravity is by 
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the use of Dirac's delta function type source with a fixed massE3). However, this 
makes Einstein's equations mathematically meaningless because of their nonlinearity. 
Physically, there is no such source in general relativity because of the existence of 
black holes. Before a body shrinks to a point, it forms a black hole whose size is 
fixed by its mass. For this reason, it has been claimed that no point particle exists 
in general relativity. 

This conclusion is not correct, however. We can shrink the body keeping the 
compactness (M/R), i.e., the strength of the internal field fixed. Namely we should 
scale the mass M just like the radius R. This can be fitted nicely into the concept of 
regular asymptotic Newtonian sequence defined above because there the mass also 
scales along the sequence of solutions. In fact, if we take the masses of two stars as M, 
and the separation between two stars as L, then ~ M/L. Thus the mass M scales 
as if we fix the separation. In the above we have assumed that the density scales 
as to guarantee this scaling for the mass understanding the size of the body fixed. 
Now we shrink the size as to keep the compactness of each component. Then 
the density should scale as e~^. We shall call such a scheme as strong point particle 
limit since the limit keeps the strength of internal gravity. The above consideration 
suggests the following initial data to define a regular asymptotic Newtonian sequence 
which describes nearly Newtonian system with strong internal gravitycj). The initial 
data are two uniformly rotating fluids with compact spatial support whose stress- 
energy tensor and size scales as and e^, respectively. We also assume that each 
of these fluid configurations would be a stationary equilibrium solution of Einstein's 
equation if other were absent. This is neccesary for the suppression of irrelevant 
internal motions of each star. Any remaining motions are the tidal effects caused 
by the other body, which will be of order smaller than the internal self-force. 
This data allows us to use the Newtonian time r = et as a natural time coordinate 
everywhere including inside the stars. 

This choice of the data leads naturally to the introduction of the body zones Ba 
and the body zone coordinates x\ defined by 

Ba = {x^ Ix'^ - ell < eR], (2-38) 
4 = e-'(^'-ei), (2-39) 

where R is some constant, and ^\{t),A = I, II are the coordinates of the origin 
of the two stars, where we have used letters with bar to distinguish the body-zone 
coordinates from their counterparts. The scaling by means that as the star 
shrinks with respect to the coordinate x^, it remains of fixed size in the body-zone 
coordinate x^. The boundary of the body zone shrinks to a point with respect to x^, 
and expands to infinity with respect to as e — > 0. This makes a clean separation 
of the body zone from the exterior geometry generated by other star. 

9^iu = {gB)f,u + {gc-B)tiu, (2-40) 

where {gB)fj.u is the contribution from the body zones, {gc-B)fj.u from elsewhere. 

Actual calculations are most easily performed in the harmonic coordinate. Take 
two stationary solutions of Einstein's equations for the perfect fluid {T^"^ (x^) , g'^ (x^)} 
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as our initial data in the body zone. As explained above every component of the 
stress-energy tensor T^*^ has the same scaling in the inertia coordinates {t,x^) 
for a rapidly rotating star. In the body zone coordinates {t,x^), these have the 
following scalings: 

TJ^ = e^T^2 ~ (2-41) 
= e-^r| ~ e-^ (2-42) 

^ij ^ ^-Ar^ij ^ ^„8^ (2-43) 

Transformed to the near zone coordinates {t,x^), x^ = ^^ + e^x'^, they take the form 

W = TV, (2-44) 
= e^tX + v\fX\ (2-45) 
= e'rj + 2e'vif'2^ + v\v\f\\ (2-46) 

where v\ = d^\/dT is the velocity of the origin of the body A. If there were only one 
body, these data would produce a stationary solution in the body zone, which moves 
with uniform velocity v\ in the near-zone coordinates. Now we know the ordering of 
the source, we can solve Eq.( p-29| ),iteratively as in the weak field case. The difference 
is that we transform the integration variables to the body zone coordinates in the 
expression of the fields whose contributions come from the body zone. 

h^j^{e;T,x^)=Ae^Y. I 'i'yA\xA - e^yA\-^A^^''{T - e\xA - e^y\,y\e), (2-47) 

A J B A 



A 



where Expanding these expressions in terms of e, we obtain 



m{e-T,x') 


= 4e^J 


- Ma 




(2-48) 


hS{e;T,x') 


= 4e^} 


Ja 




(2-49) 


~4{e;T,x'') 




ZjA 


+ Se"* 2^ -T] + 4^ II + 0(6 ), 

A A 


(2-50) 


where 














Ma 


= lime^ / d3i/yl7(e;r,y), 


(2-51) 






•Ja 


= lim€^ [ d^yAj{e;T,y), 
^^0 Jba 


(2-52) 






^A 


= lime^/ d^A%-T,v). 

^^0 Jba 


(2-53) 



The Ma defined above is the conserved ADM mass the body A would have if it 
were isolated. Without loss of generality one can set the linear momentum J\ of 
the overall internal motion of each body to zero by choosing appropriate origin of 
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the coordinates. If we did not assume the internal stationarity in the initial data, 
then Z^^ would be finite and no approximation would be possible. However under 
the condition of internal stationarity one can show that vanishes 0\ Above 
expressions for the body zone contribution are used to calculate the pseudotensor in 
the near zone and then to evaluate the contribution h^_^ outside the body zone. 
As the result we shall obtain the metric variables up to O(e^). 

r^(r, x') = 46^ ^ ^ + 0(6^), (2-54) 

h-\T, x') = 4e^ y + +2e^ V + 0{e^), (2-55) 

A l^^l A l^^l 

h^^ir, x') = 46^ V - 2e'l^2:' + 0{e'), (2-56) 

\xa\ 



A 



where 



Mj^' = lime^2 J^^d^yf4^{e;T,y), (2-57) 

Cb = T.^ACiMAJ\. (2-58) 

A 

These are the angular momentum of the body A and the quadrupole moment for 
the orbital motion, respectively. This kind of calculation is also performed at 2.5PN 
order to get the standard quadrupole formula except that the mass in the definition 
of the quadrupole moment is replaced with the ADM massEj\ The above method 
has been also applied for the calculation of spin precession and the same form of the 
equation in the case of weak gravity is obtained. Thus we can extend the strong 
equivalence principle for bodies with strong internal gravity to the generation of 
gravitational waves and the spin precession. It is an open question if the strong 
point particle limit may be taken to lead a well defined equation of motion at higher 
post-Newtonian orders. 



§3. Post-Newtonian Approximation in the (3+1) Formalism 

In the evolution of binary neutron stars, the orbital separation eventually be- 
comes comparable to the radius of the neutron star due to radiation reaction. Then, 
the point particle picture does not apply and each star of the binary begins to behave 
as a hydrodynamic object because of tidal effect. As described in the introduction we 
need not only a hydrodynamic treatment, but also general relativistic one in order to 
study the final stage of the evolution of binary neutron stars. Full general relativistic 
simulation will be a direct way to answer such a question, but it is one of the hardest 
problem in astrophysics. Although there is much progress in this directionl£3^ , it will 
take a long time until numerical relativistic calculations become reliable. One of the 
main reasons for this would be that we do not know the behavior of geometric vari- 
ables in the strong gravitational field around coalescing binary neutron stars. Owing 
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to this, we do not know what sort of gauge condition is useful and how to give an 
appropriate general relativistic initial condition for coalescing binary neutron stars. 

The other reason is a technical one: In the case of coalescing binary neutron 
stars, the wavelength is of order of A ~ vriJ^/^M"-'^/^, where R and M are the orbital 
radius and the total mass of binary, respectively. Thus, in numerical simulations, we 
need to cover a region L > A oc R^^'^ with numerical grids in order to attain good 
accuracy. This is in contrast with the case of Newtonian and/or PN simulations, 
in which we only need to cover a region X > L > R. Since the circular orbit rvf 
binary neutron stars becomes unstable at i? < lOM owing to the tidal effectsl^'H* 
and/or the strong general relativistic gravity©, we must set an initial condition of 
binary at i? > lOAf. For such a case the grid must cover a region L > A ~ lOOAf in 
numerical simulations to perform an accurate simulation. When we^siime tncoiier 
each neutron star of its radius ~ 5Af with ~ 30 homogeneous gridsE3''L3)'o)'E3''E5l', 
we need to take grids of at least ~ 500^ , but it seems impossible to take such a large 
amount of mesh points even for the present power of supercomputer. At present, 
we had better search other methods to prepare a precise initial condition for binary 
neutron stars. 

In the case of PN simulations, the situation is completely different because we 
do not have to treat gravitational waves explicitly in numerical simulations, and 
as a result, only need to cover a region at most L ~ 20 — 30M. In this case, it 
seems that ~ 200^ grid numbers are enough. Furthermore, we can take into account 
general relativistic effects with a good accuracy: In the case of coalescing binary 
neutron stars, the error will be at most ~ M/R ~ a few x 10% for the first PN 
approximation, and ~ (M/R)'^ ~ several % for 2PN approximation. Hence, if we 
can take into account up through 2PN terms, we will be able to give a highly accurate 
initial condition (the error < several %). 

For these reasons, we present in this section the 2.5PN hydrodynamic equations 
including the 2.5PN radiation reaction potential in such a way that one can apply the 
formulation directly in numericaJLsimulations. As for the PN hydrodynamic equation, 
Blanchet, Damour and SchaferllB* have already obtained the (1+2.5) PN equations. 
In their formulation, the source terms of all Poisson equations take nonvanishing 
values only on the matter, like in the Newtonian hydrodynamics. Although their 
formulation is verv useful for PN hydrodynamic simulations including the radiation 
reaction they did not take into account 2PN terms. In their 

formulation, they also fixed the gauge conditions to the ADM gauge, but in numerical 
relativity, it has not been known yet what sort of gauge condition is suitable for 
simulation of the coalescing binary neutron stars and estimation of gravitational 
waves from them. First, we develop the formalism for the hydrodynamics using the 
PN approximation. In particular, we use the (3+1) formalisin_pf general relativity 
so that we can adopt more general class of slice conditions CJ^'cP. Next, we present 
methods to obtain numerically terms at the 2PN order 

3.1. (3+1) Formalism for the post- Newtonian Approximation 

According to the above discussion, we shall apply IheJ 3+1) formalism to formu- 
late the PN approximation. In the (3+1) formalism EP'EHP'Ea^ spacetime is foliated 
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by a family of spacelike 3D hypersurfaces whose normal one-form is taken as 

= i-a, 0). (3-1) 

Then the line element takes the following form 

ds^ = - (a^ - PiP')dt^ + 2Pidtdx' + jijdx'dx^ , (3-2) 

where a, and jij are the lapse function, shift vector and metric on the 3D hyper- 
surface, respectively. Using the (3+1) formalism, the Einstein's equation 

Gnu = SttT^,,, (3-3) 

is also split into the constraint equations and the evolution equations. The formers 
are the so-called Hamiltonian and momentum constraints which respectively become 

tri? - KijK'^ + K^ = IGirpH, (3-4) 

DiK'j-DjK = 8TrJj, (3-5) 

where Kij, K, iiR and Di are the extrinsic curvature, the trace part of Kij, the 
scalar curvature of 3D hypersurface and the covariant derivative with respect of jij. 
PH and Jj are defined as 

PH = T^uu^h", (3-6) 

The Evolution equations for the spatial metric and the extrinsic curvature are re- 
spectively 

—-fij = -2aKij + DiPj + Djfii, (3-8) 



+rDmKij - 87ra [s^j + ^7i^- {pn - S\)^ , (3-9) 
^7 = 27(-ai^ + A/30, (3-10) 



-K = a{trR + K^) - D'Dia + ^WjK + 47ra(5', - 3ph), (3-11) 



d_ 

w 

where Rij, 7 and Sij are, respectively, the Ricci tensor with respect of 7ij, the 
determinant of jij and 

Sij = na\l'r (3-12) 

Hereafter we use the conformal factor -0 = 7^/^^ instead of 7 for simplicity. 

To distinguish the wave part from the non-wave part (for example, Newtonian 
potential) in the metric, we use 7^^ = ■i/'~^7ij instead of 7ij. Then det(7ij) = 1 is 
satisfied. We also define as 



Aij = ^-^Aij = [Xij - \li3K^ . (3- 13) 
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We should note that in our notation, indices of Aij are raised and lowered by 7ij, 
so that the relations, A^j = A^j and A''^ = ip^A^^ , hold. Using these variables, the 
evolution equations (3-8)-(3-ll) can be rewritten as follows; 



dn 



~ df3^ df3^ 



-A- 

dn 



1 

2 



2, dp^ 



(3-14) 



2 a/?" 



+a{KAij - 2AiiA + -Q—^mj + -Q^^rai - ^-^^iJ 



dn 
d 



-aK + 



1 



ih 



1 



(3-15) 
(3-16) 



i-K = a{ AijA'^ + ) - -T^" - 75 7^V,ika,i + 47ra(5\ + pn), (3-17) 



dn" ' 3" J ilA 

where Di and A are the covariant derivative and Laplacian with respect to and 

(3-18) 

The Hamiltonian constraint equation then takes the following form. 



— = — -If — 
dn dt dx^ 



All) = ^tri?V - 27rpffV^ - ^ ( Aj^'^ - ] , 



(3-19) 



where tri? is the scalar curvature with respect to 7ij. The momentum constraint is 
also written as 



(3-20) 



Now let us consider Rij in Eq.(3-15), which is one of the main source terms of 
the evolution equation for Aij. First we split Rij into two parts as 



Rij — Rij + R!^j , 



(3-21) 



where Rij is the Ricci tensor with respect to 7ij and i?^- is defined as 
R% = —DiDj^ - ^lijD'^Dki' + ^{Dii^){Dji^) - ^%{Dk^){DH)- (3-22) 

Using det(7jj) = 1, Rij is written as 

1 
2 



R^ 



1] 



hij^kk + hji^ii + hii^ij + f^\j.{hij^i + hii^j — hij^i) 



+f'^\hkj,ii + hkiji — hij^_ 



kl) 



T^k 



(3-23) 
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where denotes d/dx^, F^^ is the Christoffel symbol, and we spht and 7*-' as 
^ij + hij and (5*-^ + /^-^ , by writing the flat metric as 5ij . 

We shall consider only the linear order in hij and fij assuming \hij\, \fij\ <^ 1. 
(As a result, hij = —P^ + 0{h?').) Such a treatment is justified because hij turns 
out to be 2PN quantity in our choice of gauge condition (see below). Here, to clarify 
the wave property of 7jj, we impose a kind of transverse gauge to hij as 

hij,j = 0, (3-24) 

which was first proposed by Nakamura in relation to numerical relativity il^. Here- 
after, we call this condition merely the transverse gauge. The equation (3-14) shows 
that this condition is guaranteed if /?* satisfies 

- = (-2aiii + + - \%f3'^ . (3-25) 

Using the transverse gauge, Eq.(3-23) becomes 

Rij = -]^Ahij + 0{h'^), (3-26) 

where A is the Laplacian with respect to 5ij. Note that tri? = 0{h'^) is guaranteed in 
the transverse gauge because the traceless property of hij holds in the linear order. 

We show the equations for the isentropic perfect fluid ( ^-161 ). The conservation 
law of mass density, (|^), may be written as 

where p^, is the conserved density defined as 

= aip^pu^. (3-28) 

The equation of motion and the energy equation are obtained from the conservation 
law ( |2-18| ) which takes the following forms. 

^ + = -o^^'P. - c^o^.S' + S,P\, - ^,S,S,^',, (3-29) 



and 



where 



Si = aip^ip + pn + p)u'^u, = pji + n + - )ui{= iij^Ji), 



= a^^{p + pn + P){u' 
H = aip^pHu^ = p^n, 



0\2 



P 

{PH + pW 



a 



-^-^ = -/3' + ^- (3-31) 

Finally, we note that in the above equations, only /?* appears, and (3i does not, so 
that, in the subsequent section, we only consider the PN expansion of not of (5i. 
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3.2. Post- Newtonian approximation in the (3+1) formalism 

Next, we consider the PN approximation of the above set of equations. First 
of all, we review the PN expansion of the variables. Each metric variable may be 
expanded up to the relevant order as 

^ = 1 + (2)^- + {4)V' + {6)V' + {7)V' + • • • , 
a = 1 + (2)0 + (4)0 + (6)0 + (7)a + . . . , 

1-U+' 



+ X] + (6)0 + (7)0 + . . 



hi 



A 



K 



(3-32) 



(4) hij + (5) hij + . . . , 

(3)ir+(5)i^ + (6)K + ..., 

where subscripts denote the PN order(e"') and U is the Newtonian potential satisfying 

AU = -iirp. (3-33) 

X depends on the slice condition, and in the standard PN gauge0\ we usually use 
^ = —X/2, which satisfies 



/ 1 3 P 

A'P = -ATrpi + U + -n + -- 
V 2 2 p 



(3-34) 



Note that the terms relevant to the radiation reaction appear in (7)0, (s)/?* and 
[5)hij, and the quadrupole formula is derived from (7)0 and (j,)hij. 
The four velocity is expanded as 



u 



Ut 



1 



1_ 



l + e'i-v' + U +e''l-v^ + -v'U + -U' + ^,)PV-X +Oie 



u = V 



l + e'{ -v' - U ] + -v"" + -v'U + -U' + X ] + 0{e 



+ 0{e') 



W ilv^ + Iv^U + AU^-X + 4(4)^ + (3)/3-''r;^ 



+ e ( (6)/3 + {5)/ij, 



(3-35) 



where = 0(e) and = v'^v'^ . From u^u^j^ = —1, we obtain the useful relation 



{au^f = 1 + Y^UiUj 



1 + + (f^ + 4v^U + 2(3)/3*?;^) + O(e^). 



(3-36) 
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Thus ph, Ji and Sij are respectively expanded as 



PH = e p 
Ji = e'p 
Sij = e^p 



I + v"^ + n ] + e*{v* + \ AU + n + 



4 J „,4 



P 



vHl + e^ [v^ + 2,U + n + 



P 
P 



P 



UP , 
P 



P 



S, 



+ 0{e% 



+ 0(e^). (3-37) 



The conformal factor ip (and a in the conformal slice) is determined by the 
Hamiltonian constraint. In the PN approximation, the Laplacian with respect to 7*-' 
for the scalar is expanded as 



A = A- (e^(4)/i,, + e\5^hij)d,dj + 0{e^) 
At the lowest order, the Hamiltonian constraint becomes 

^{2)V' = -27rp. 



(3-38) 



(3-39) 



Thus, (2)CK = — 2(2)V' = ~U is satisfied in this paper. At the 2PN and 3PN orders, 
the Hamiltonian constraint equation becomes, respectively. 



Z\(4)V' = -2ttp[v^ + n + -u 



(3-40) 



and 



Z^(6)^ = -27Tpi^V^ + v-'(^n+^ + y + 2(3)/?^^;^ + + + 5(4) 

(3-41) 



+i^{4)hijU^ij 



{3)^ij{3)^ij 



[(3) 



The term relevant to the radiation reaction first appears in ^^^ip whose equation 
becomes 

1 



^(7)^ 



-,{5)hijU^ij. 



(3-42) 



Hence, (7)a may be also relevant to the radiation reaction, depending on the slice 
condition. 



Let us now derive equations for p\ From Eq.( p-25D , the relation between {3)Aij 
and (3)/3* becomes 



(3-43) 



where we used the boundary condition that (3)^ij and (3)/?* vanishes at the spatial 
infinity. (3)^jj must also satisfy the momentum constraint. Since (3)^ij does not 
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contain the transverse-traceless (TT) part and only contains the longitudinal part, 
it can be written as 



(3)^ji = (3)^i,i + (3)^i,j - -^^ij(3)Wk,k , 



(3-44) 



where (3) Wj is a vector on the 3D hypersurface and satisfies the momentum constraint 
at the first PN order as follows; 



1 2 

^{3)Wi + -(3)Wjji - -(3)K^i = 8TTpv\ 

From Eq.( 3-43| ), the relation, 

(3)/3' = 2(3) TF, , 
holds and at the first PN order, Eq.(3-16) becomes 

3f7 = -(3)/^ + (3)/?',/ , 



(3-45) 
(3-46) 
(3-47) 



where U denotes the derivative of U with respect to time. Thus Eq.( 3-45 ) is rewritten 



as 



(3-48) 



This is the equation for the vector potential at the first PN order. 

At the higher order, is also determined by the gauge condition, hijj = 0. 
The equation for (5)/3* is obtained by using the momentum constraint and the 2PN 
order of Eq.(3T6) as follows. 



2?7 + 77 + - ) + (3)/3* 
1 



8f/,,(3)A 



+(5)7C,i - C/(3)7r,, + -U^i^)K - 2(4)^,, + -(UU),, + ((3)/3't/,0,i • 

(3-49) 

The equation for (6)/3* takes the purely geometrical form since the material contri- 
bution Jj at the 1.5PN order vanishes. 

^(6)/9' = (6)K,^ . (3-50) 

Then, let us consider the wave equation for hij. From Eqs.(3-14), (3-15), (|^) 
and (3-26), the wave equation for hij is written as 

nhii = ( 1 - — ]Ahij + { — - — ]h 



2a 



DiDj - -^ijA ) a - ^ ( D.i^Dja + Dj^jDia - -^D^i^DkU 



-167r- 



a 



4 \ 



dn 
(3-51) 
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where □ is the flat spacetime wave operator defined as 

92 



□ 



dt^ 



+ A. 



We should note that (4)Tij has the TT property, i.e., (4)T, 
is a natural consequence of the transverse gauge, h. 



(3-52) 



and (4)Tij = 0. This 



and ha = 0{h^). Thus 



(4)hij is determined from 



A.(A)hij = {4)Tij- (3-53) 

Dmce 0{h^) turns out to be O(e^), it is enough to consider only the linear order of 
hij in the case when we perform the PN approximation up to the 3.5PN order. We 
can obtain (5)/iij by evaluating 



]_d_ 

An dr 



(3-54) 



and the quadrupole mode of gravitational waves in the wave zone is written as 



hlf{T,^) = -^ lim 



1 



4tt |> 



{4)Tii(r-e|x-y|,y) 

— j j ay. 



(3-55) 



Later, we shall derive the quadrupole radiation-reaction metric in the near zone using 
Eq.(|3^. 

Finally, we show the evolution equation for K. Since we adopt slice conditions 
which do not satisfy K = (i.e. the maximal slice condition), the evolution equation 
for K is necessary. The evolution equations appear at the IPN, 2PN and 2.5PN 
orders which become respectively 



d / n P 

— (3)/C = ATTpi^2v^ + n + 2U + 3- 



dr 
d 

— (5)i^ = Airp 



AX, 



2t;^ + v^ 6C/ + 277 + 2- 



P 



+4(3)/?' 



+ {3)Aj{3)^ij + 3(3) 



n + 



3P 



(3-56) 



P J 

(4)hijU^ij + (3)/?* (3)7^,1 



-A(7)a 



(3-57) 
(3-58) 



We note that for the PN equations of motion up to the 2.5PN order, we need 

(2)0, (4)0, (6)«, (?)«> (2)V') (4)'0) {3)P'\ (5)/^*, (4)^ji, (5)^ii, (5)K and (6)K. 

Therefore, if we solve the above set of the equations, we can obtain the 2.5 PN 
equations of motion. Up to the 2.5PN order, the hydrodynamic equations become 

+p, \uji + n + - + h^ -u + lv^ + iv'^u 

L I P 2 8 
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+ {-v' -u]{n + 



p 



P 

-x,Ai + n + - + - 
p 2 



+ '^V^(4.)tp,i - (6)«,j - (7) 



P V 



+^;^ (3)/3^ i 1 + 77 + - + - + 3C/ + (5)/?^ i + (g)/?;' 



P 2 



+ 0(6^), (3-59) 



dH d{Hv^) 
dr dx^ 



-P 



dt\2 



v'^ + 3U)v^ } + 0{e 



(3-60) 



where we have omitted e" coefficients for the sake of simphcity, and used relations 



Si = p* 



P 



P 



i + n + - + — + — [n + j]+^v^ + 2v'^u + (3)/5-'f-'^ 



p 2 2 
l + i7 + ^ + ^+3C/) +(3)/3^ 



+ 0(6'^), 

(3-61) 



3.3. Strategy to obtain 2PN tensor potential 

Although the 2PN tensor potential is formally solved as 



{A)hij{T, x) 



1 

47r 



(4)^ii(T, y) 3 



|x - y| 



(3-62) 



it is difficult to estimate this integral numerically since (4)Tjj — > 0(r~ ) for r — > oo 
and the integral is taken all over the space. Thus it is desirable to replace this 
equation by some tractable forms in numerical evaluation. We shall present two 
methods to do so. One is to change Eq.( |3-62j ) into the form in which the integration 
is performed only over the matter distribution like as in the Newtonian potential. 
The other is to solve Eq.(3-62) as the boundary value problemo'. 

Strategy 1: Direct integration method 

The explicit form of (4)Tjj is 



(4)rij = -2dij{X + 2(4)V^) + UdijU - 3U,iUj + 6ijU,kU,k - 167r( 

2 



where 



pv'v^ - -Sijpv'^ 



52 1 

dij = ^ ■ „ — r 6ijA. 

ox'^ox^ 3 



(3-63) 



(3-64) 



Although (4)Tjj looks as if it depends on the slice condition, the independence is 
shown as follows. Eq.( |3-4^ ) is solved formally as 



1 

47r 



(3)^ ^3 



|x-y| 



d'y 



(3-65) 
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where 



Pi 



pv 



|x-y| - 2 
From Eqs.( |340D and (3-57), we obtain 



P U 



(3)i^ = -A{X + 2(4) V') + 47rp[v^ + 3 



Combining Eq.( |3-65| ) with Eg. ( 3-67 ), the equation for (3')/3* is written as 



(3-66) 



(3-67) 



i3)P' =Pi-iX + 2(4)V),i 



(^pv^ + 3P- pU/2^ 
J |x-y| 



Using this relation, the source term, (4)rjj, is spht into five parts 



-iS) 



where we introduced the fohowing notations. 

1 

3 



(V) 

ij ' 



(3-68) 



(3-69) 



iS) 



-levrf 



pv'-v^ - ^Sijpv'^ ] , 



1 



(4)T, 



(p) 



(4)T, 



/ I 1 ^ y + / I r ^ 2/ - o^'. 

oxJ J |x — y| ax* j |x — y| 3 

/ /j|x-y|(i^y, 

pv"^ + 3P - pC//2 



,k\- 



-d y, 



|x - y| 



2a 



|x - y| 



d'y. 



(3-70) 



Thus it becomes clear that (4)/ijj and (5)/iij as well as (4)Tjj are expressed in terms 
of matter variables only and thus their forms do not depend on slicing conditions, 
though values of matter variables depend on gauge conditions. 

Then, we define (4)/iSf = /\"\4)Tif\ {A)h[J^ = Z\-\4)4^\ {A)h\f^ = ^ 



(4)% - ^ (4)Tij and (4)ft,jj. 



{4)T. 



(C) 

ij ! 



A ^(4)T^j^\ and consider each term separately. 



First, since (4)Tjj is a compact source, we immediately obtain 

,,2 



pv'v^ - \5ijpv' 



|x - y| 



Second, we consider the following equation 

Z\G(x,yi,y2) 



|x-yi||x-y2| 



(3-71) 



(3-72) 
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It is possible to write (i)hf^^ using integrals over the matter if this function, G, is 
used. Eq. (|372|) has solutions i^'O, 

G'(x,yi,y2) = ln(ri + r2 ± ri2), 

where 



(3-73) 



n = |x - yil, r2 = |x - y2|, ri2 = |yi - y2| 



(3-74) 



Note that ln(ri + r2 — ri2) is not regular on the interval between yi and y2, while 
ln(ri + r2 + ri2) is regular on the matter. Thus we use ln(ri + r2 + ri2) as a kernel. 
Using this kernel, UU^ij and U^iUj are rewritten as 



UU, 



dx^dx^ \J |x 



^^^^^ r^^^i 



yil 



^^^^^ r'i^^. 



|x - y2| 



a2 



1 



9yfay^ V|x-yi||x-y2| 
92 



A / d^yi(i^y2p(yi)p(y2) ^ , ^ + r2 + ru), 

•1 ov\ ovi 



^ f ^(^l) d'y. 



dx 



|x - yil 



^^^^^ 



|x-y2| 



d^yid^y2p{yi)p{y2) 



dy\dy{ 

(— 

dy{dy{ V|x - yi||x - y2| 



A / (i^yi(i'y2P(yi)p(y2)T— TT Hn + ^2 + ri2). (3-75) 
J oy\dy^2 



Thus we can express (4)/ij-j using the integral over the matter as 



{4)%- 



(U) 



d yid y2p{yi)p{y2] 



( 



\dy\dy{ 3 
X ln(ri + r2 + ri2), 
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dy\dyl 3 



'5ijAi2 



(3-76) 



where we introduced 



Ai 



92 



QykQyk 



A 



92 



12 



(9^2 



(3-77) 



Using relations Z\|x — y| = 2/|x — y| and Z\|x — yl'^ = 12|x — y|, (^4-^h'l^\ (4)/i|j^ and 



h]- are solved as 



2^ j{pv^y\^-y\d^y + 2-^ J (pv'y\^ - y\d^y + ^5i, J p\^ - y\d%, 

(3-78) 



(4) 



(P) ^ 1 

'J 12 dx'dx^ 



jjp\^- yfct'y - l^ii J P|x - y\d^y, (3-79) 
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(4)' 



pv'^ + 3P- pU/2 



X ■ 



(3-80) 



Thus we find that the 2PN tensor potentials can be expressed as the integrals only 
over the matter. 

Strategy2: Partial use of boundary value approach 

Although the above expression for (4)/iij is quite interesting and might play 
an important role in some theoretical applications, it will take a very long time 
to evaluate double integration numerically. Therefore, we propose another strategy 
where Eq.( |3-53| ) is solved as the boundary value problem. Here, we would like to 
emphasize that the boundary condition should be imposed at r(= |x|) 3> |yi|, |y2|, 
but r does not have to be greater than A, where A is a typical wave length of 
gravitational waves. We only need to impose r > i? (a typical size of matter). This 
means that we do not need a large amount of grid numbers compared with the case 
of fully general relativistic simulations, in which we require r > X ^ R. 

First of all, we consider the equation 



A 



(S) , (U) 



(3-81) 



Since the source (4)7"/^'^^ behaves as 0(r ^) at r ^ oo, this equation can be accurately 
solved under the boundary condition at r > i? as 



iS) 



-I, 



where 



+0{r' 



^ijk 



(3-82) 



px^x^x^d^x, 
Sijk = I piv'-x^ - v^x'')x'^dPx. 



(3-83) 

Next, we consider the equations for {^^^h^^\ {A)h\j^ and (4)/i-^^ Using the iden- 



tity, 

p = -{pv% = (pvV),,, + AP- (pUi),, 
we find the following relations; 



(3-84) 



p\x - y\d y 



p|x-y|3d=^y = 3 / d-'y 



x-y| 

■ (x* -y'){x^ - y^) 



pv V 



X ■ 



+ [AP + pv^ - pU^i{x' -y')]\^-y 



(3-85) 
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Using Eqs.(3-85), (4)h^^\ (4)h\f and ^^^h^j > in Eqs.(|37^)-(|3^ can be written as 



iv) 



|x-y| 



|x - y| 



|x-y| 

(3-86) 



A: A; 



(4) ii 4 9a,ia3;i 



X — y| ox^ J |x — y| 



2 



k d 



|x-y 



|x-y| 



+ . (3.87) 



dx3 



|x - y| 



and 



pv'^ + 3P- 



|x- y| 



(3-88) 



where P' = P + pv'^ /A + pU^ii/ /A. Furthermore they are rewritten in terms of the 
potentials as follows. 



(4)/.Sf = 2{x\,)Pi + X^(3)P^ - g.,) + ^5., - X\,^P'^ ) , 



- 1 

1 r 5 



Pk 



Q 



kk 



1 



5x* 



V,. 



ki J 



(4)% 



^{Qfx^ + Qfx 



Qi'J-Qfj) + lQ^'^^^^^ (3-89) 
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where the potentials are defined as 

AQij = -iTr(x^{pv'y -\-x\pv^ 



AQ 



(I) 



AV, 



ifw) 

't 



AV^ 
AV, 



(P) __ 
(pu) 



AV, 



(pU) 



-Att(^PV^ + 2,P-]^pU^, 

— 47r/>w*f-', 

—Airpy'^v-'x^ 

—47rp{v^x^)'^, 

-AttP', 

-AttP'x\ 

—AnpU^iX^ . 



(3-90) 



It should be noted that these Poisson equations have compact sources. 

In this strategy {A)h[j^ and {4)hf^^ are solved as the boundary value problem, 
while other parts are obtained by the same method as in Newtonian gravity. 

Strategy 3: Boundary Value approach 

Instead of the above procedure, we may solve the Poisson equation for (4)/iij as 
a whole carefully imposing the boundary condition for r ^ i? as 



-'fcfe g ^ki ^ ^"v^kk g"u" " ^ki r 



+- 



2 . 4 



+ [-^n iS,k,+S,u)-^{n^S,kk + n^S,kk) 

2 

+2n''n\n'Sjki + Siki) + 2n'n^n^Skii + -SijU^Sku 



+ 0(r-3). 

(3-91) 



It is verified that 0{r~^) and 0{r~'^) parts satisfy the traceless and divergence-free 
conditions respectively. It should be noted that (4)/iij obtained in this way becomes 
meaningless at the far zone because Eg. ( 3-53 ), from which (4)/ijj is derived, is valid 
only in the near zone. 
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3.4. The Radiation Reaction due to Quadrupole Radiation 

This topic has been already investigated by using some gauge conditions in 
previous papers . However, if we use the combination of the conformal shce 

and the transverse gauge, calculations are simplified. 

3.4.1. conformal slice 

In combination of the conformal sliced and the transverse gauge, Eq.( p-5^ ) 
becomescP 



-167r( pv^v^ — -6ijpv^ 



+ uu 



o ^3 



(3-92) 



lines becomes —2f\^^ and the third line becomes 6rl-\^^ /5, where -J-^^^ = d^Fij/dt^ 
and 



From a straightforward calculation, we findo' that the sum of the first and second 

1 

-I] — ^ij — -L 

Thus, (5)/iij in the near zone becomes 



Fi 



I, 



(3-93) 



1 r(3) 
5 ' 



(3-94) 



Since hij has the transverse and traceless property, it is likely that (5)/ijj remains 
the same for other slices. However it is not clear whether the TT property of hij 
is satisfied even after the PN expansion is taken in the near zone and, as a result, 
whether (^)hij is independent of slicing conditions or not. The fact that slicing 
conditions never affect (5)/ijj is understood on the ground that (4)rjj does not depend 
on slices, as already shown in Eq.(3-69). 

Then the Hamiltonian constraint at the 2.5PN order, Eq.( 3-42| ), turns out to be 



5 



where x is the superpotential0^ defined as 

j p\y^-y\d^y 



X 



(3-95) 



(3-96) 



which satisfies the relation Ax 
form, 



-2U . From this, we find (j-^'ip takes the following 



(7)V' 



5 

5 



^d'y 



(3-97) 
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Therefore, the lapse function at the 2.5PN order, (y^a = — 2(7)^/^, is derived from U 
and Ur, where Ur satisfies S 

AUr = (3-98) 

Since the right-hand side of Eq.(3-58) cancels out, (e)-^ disappears if the (e)-?^ does 
not exist on the initial hypersurface, which seems reasonable under the condition that 
there are no initial gravitational waves. Also, (6)/9* vanishes according to Eq.( |3-5Ci| ). 
Hence, the quadrupole radiationreaction metric has the same form as that derived 
in the case of the maximal slice 

From Eq.( |3-29| ), the PN equation of motion becomes 

yi + E± + U i + i^^^^ + Ff""^ + F^-^""^ + ©(c^), (3-99) 

p 

where F^^^ and Ff^^ are, respectively, the IPN and 2PN forces. Since the radi- 
ation reaction potentials, {j,)hij and (7)Q;, are the same as those by Schafer (1985) 
and Blanchet, Damour and Schafer (1990) in which they use the ADM gauge, the 
radiation reaction force per unit mass, F^-^p^ = , is the same as their force and 



Fl = - {{{->)hijV^y + v^v^f^(^5)hij + (7) 



^rff^y + g^S'^v, + g^r^?; ''('•^) ' |.-;|3 " ^'y, 

(3-100) 

Since the work done by the force (3-100) is given by 
W = j pF[v'd^x 

we obtain the so-called quadrupole formula of the energy loss by averaging Eq. (3-101) 
with respect to time as 



(^) = -^(4^)^(^)-) + 0(e«). (3-102) 



3.4.2. Radiation reaction in other slice conditions 

In this subsection, we do not specify the slice condition. In this case, the reaction 
force takes the following form 

F[ = ^{Fifv^- + lAfv% - (7)a i - (6)/?^ + v\,)P^ , - v\,)P^ (3-103) 



Here, (7)0 corresponds to the slice condition. From Eq. ( 3-103 ), we obtain the work 
done by the reaction force as 



W = j pF[v'd^x 
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d^xpv^ 



(3-104) 



Explicit calculations are done separately: For the first and second terms of 
Eq.(3-104), we obtain 



4 d 
5di 



4- 



(3) 



d'^xpv^v-' 



_^ 7:(3) 
5 

5 



d^xp{x)v 
d^xpv^v^ , 



1 



(3) r(3) 



5 *J «J 



d 
dx^ 



,3 p{y){x' - y'){x^ -y^) 



and 



d^xp^lfv^iv^v' 



2 r{3) 
5 



d^xpv^vK 



As for the fourth term in the integral of Eq.(3-104), we find the relation 



(6)/?' 



d 



[x^ — y^){x^ — y^) 



X 



(3-105) 
(3-106) 

(3-107) 



which is given by Eqs.( |3^ , (3- 58) an d ( ^9^)- 

Using Eqs.(3-104), (3-105), (l3a06| ) and (|3-107| ), we obtain 



-(3) 



pv^v^d'^x 



i 7:(3)_7:(3)ii 
5 



(3-108) 



This expression for W does not depend on the slice condition. However, this never 
means that the value of W is invariant under the change of the slice condition, since 
the meaning of the time derivative depends on the slice condition. 

3.5. Conserved quantities 

The conserved quantities are gauge- invariant so that, in general relativity, they 
play important roles to characterize various systems described in different gauge 
conditions. From the practical point of view, these are also useful for checking the 
numerical accuracy in simulations. Thus, we show several conserved quantities in 
the 2PN approximation. 

Conserved Mass And Energy 

In general relativity, we have the following conserved mass; 



p^d^x. 



(3-109) 



In the PN approximation, defined by Eq.( 3-28| ) is expanded as& 
1 



+ [^v^ + ^v'U + ^U' + 6(4) V + i3)(3'v' ] + (6)5* 



+ 0(e^), (3-110) 



Post-Newtonian Approximation 



31 



where (6)5* denotes the 3PN contribution to p*. This term (g)*^* will be calculated 
later. 

Next, we consider the ADM mass which is conserved. Since the asymptotic 
behavior of the conformal factor becomes 



V' = 1 + 



Madm 
2r 



(3-111) 



the ADM mass in the PN approximation becomes 



Madm 



1 

5 



AipcFx 



13 



P 



1+ [v' + n + -u] + [v^ + —v'u + v'n + -v^ + -un 



+-C/^ + 5(4)V + 2(3)/3V + 



+{&)5adm + 0{e' 



(3-112) 



where (^q-^Sadm denotes the 3PN contribution. This term (6)5a_dm will be calculated 
later. 

Using these two conserved quantities, we can define the conserved energy as 
follows; 



E = Madm - M* 



d^xp 



+ 



1 



-v^ + n - -u 



1 



IGvrp 

En + EiPN + E2PN + 



(3)^ji{3)^ji - ^{3)K ) + ( (g)Sadm - (6)^*] + 0{c ) 



(3-113) 



We should note that the following equation holds 



(3)^1^,(3) Ai.d^x = -8^ J pv\^)l3'd'x + J [^i3)K' + 2U^s)Kjd'x, (3-114) 

where we used Eqs.( p-47| ) and ( |3-48| ). Then, we obtain the Newtonian and the first 
PN energies as 



En = I pl^^v^ + n-^uy^ 



X, 



(3-115) 



and 



ElPN 



d^x 



P + -v'u + v'n + -v' + 2un - -u' + -(3)/3* 



(3-116) 
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EiPN can be rewritten immediately in the following form used by Chandrasekhar 



23L 



ElPN 



d'^xp 



-v'^ + -v^U + [ n + 



P 



2 2 ^ 



(3-117) 



where qi is the first PN shift vector in the standard PN gaugec^), which turns out 
to be (3)-?^ = in the (3+1) formalism and it satisfies 



(3-118) 



The total energy at the 2PN order E2PN is calculated from the 3PN quantities 
(6)5* and (1^)5 ADM- In a straightforward manner, we obtain 



(3-119) 



where 



(6)' 



^ ,6 , ^^„,47 



16 



93 



+ -V'^U + 5(4)V + -^U' + -^3)P'V' -X + 6(6)V + 15C/(4) V 



Next, we consider (^q)6adm- The Hamiltonian constraint at O(e^) becomes 



(3-120) 



32 



'2(4)/lfc«,m(4)^fcm,/ + {4)hkl,m{A)hkl,rri) — '^T^{6)P4' 



16 



(3)^ij{3)^ij - 3(3)^ ), (3-121) 



where we define (6)Pv 

,,6 I „,4 



(6) Pip = P 



5 



{^^{n + ^) + 9(4)^ -2X + 20[/2| 



+77 ( 5(4) V + ^U'\+ 5(6) V + 10t/(4)V' + ^U' + (4)/i»jt'*t;^ 



P 13, 



+2^S)P'v'^2v' + n+- + y C/| + 2(5)/?*^* + (3)/?^3)/3^ 

Making use of the relations (4)hijj = and [e)hijj = 0, we obtain 



(3-122) 



{q)Madm = I d X(6)P^ + — / d y[ (3) ^^^(5)^^^ - -(3) 7^(5) 



+- 



327r 



d yU[ (3)Aij(3)Aij 



[(3) 



(3-123) 



where we assumed (6)^ij — > 0(?' ^) as r — > 00. Although this assumption must 
be verified by performing the 3PN expansions which have not been done here, it 



Post-Newtonian Approximation 



33 



seems reasonable in the asymptotically flat spacetime. From (g)Madm and (q)M^,, 
we obtain the conserved energy at the 2PN order 



E2PN = {6)MaDM — (6)^* 

d^xp 



+l(4)hijvV + 2(3)/3^i;^| (^77 + ^) + 5C/| + ^,)P'v' + ^(3)/3\3)/3' 

+^ / ((4)/iijC/C/,ii + {3)Aiji5)Aij - ^(3)i^(5)i^) . (3-124) 

Here we used Eq.(3-41) and the relation, / d^xp^Q-^ip = — ^ / d^xUAf^Q-^tp, in order 
to eliminate (gjV- 

Conserved linear momentum 

When we use the center of mass system as usual, the linear momentum of the 
system should vanish. However, it may arise from numerical errors in numerical 
calculations. Since it is useful to check the numerical accuracy, we mention the 
linear momentum derived from 

Pi = — lim ([( KiiU^ - Kn^dS 



Svr »'^oo 

= — lim /fv^i^n^' - -Kn^dS, (3-125) 
Svr r~*oo J \ 3 / 

where the surface integrals are taken over a sphere of constant r. Since the asymp- 
totic behavior of Aij is determined by 

(3)4- = ^((3)/?j,- + (3)/3fi - lSij{3)l3l^ + 0{r-% (3-126) 

and 

(5)4- = I ((5)/?:^ + (5)/3f. - ^5..(5)/9,z) + + 0{r-^)^ (3-127) 

the leading term of the shift vector is necessary. Using the asymptotic behavior 
derived from Eq.( 3-4S ) 

7 L 1 n'uH 



we obtain 



(3)/5^ = -o - - + 0{r-'), (3-128) 

\ / 2 r 2 r 



J ((3)/?^,- + (3)/?fi - l^^H3)(3l^^'dS = 16nh, (3-129) 
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where we defined h = J pv'^d^x. Tlierefore the Newtonian hnear momentum is 

Pn' = J d^xpv\ 
Similarly the first PN linear momentum is obtained as 

PiPN ' = j d^xp 
We obtain P2P n * by the similar procedure as 



(3-130) 
(3-131) 



P2PN ' 



d'^xpv' 



2(3) /3V + 10(4) + + 

+—u^ + muv'^ + v^-x 
4 



n + ^ 

P 



(3-132) 



§4. Formulation for Nonaxisymmetric Uniformly Rotating Equilibrium 

Configurat ions 

We now consider the construction of the spacetime involving a close binary neu- 
tron stars as an application of PN approximation. It is assumed that the binary 
stars are regarded as uniformly rotating equilibrium configurations. Here, we men- 
tion the importance of this investigation. To interpret the implication of the signal of 
gravitational waves, we need to understand the theoretical mechanism of merging in 
detail. When the orbital separation of binary neutron stars is < lOGM/c^, where M 
is the total mass of binary neutron stars, they move approximately in circular orbits 
because the timescale of the energy loss due to gravitational radiation is much longer 
than the orbital period. However, when the orbital separation becomes 6 — lOGM/c^, 
the circular orbit cannot bemaintained because of instabilities due to the GR grav- 
ityil) or the tidal fieldil^'El^. As a result of such instabilities, the circular orbit of 
binary neutron stars changes into the plunging orbit to merge. Gravitational waves 
emitted at this transition may present us important information about the structure 
of NS's since the location where the instability occurs will depend sensitively on the 
equation of state (EOS) of NSE^^'H^'llHs. Thus, it is very important to investigate 
the location of the innermost stable circular orbit (ISCO) of binary neutron stars. 

In order to search the ISCO, we can take the following procedure: First, ne- 
glecting the evolution due to gravitational radiation, we construct equilibrium con- 
figurations. Next, we take into account the radiation reaction as a correction to the 
equilibrium configurations. The ISCO is the orbit, where the dynamical instability 
for the equilibrium configurations occurs. Hence we shall present a formalism to 
obtain equilibrium configurations of uniformly rotating fluid in 2PN order as a first 
stepB*, though in reality due to the conservation of the circulation for gravitational 
radiation reaction, tidally locked configuration such as uniformly rotating fluid is not 
a good approximation. 

4.1. Formulation 

We shall use in this section the maximal slice condition iC-* = 0. As the spatial 
gauge condition, we adopt the transverse gauge j- = in order to remove the gauge 
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modes from ^ij. In this case, up to the 2PN approximation, each metric variable is 
expanded asu^ 



V' = i + 3T + :z(4)^ + 0(O, (4-1) 



2 c4 



a = l-^[/ + ^(^ + X) +-^(6)a + 0(c-^), (4-2) 

P' = ^3i3)P' + ^i,)(3' + 0{c-'), (4-3) 
1 



= ^ij + ^/lii + 0(c (4-4) 

In this section, we use 1/c instead of e as the expansion parameter because of its 
convenience in numerical applications. 

For simplicity, we assume that the matter obeys the polytropic equation of state 
(EOS); 

p = {r- i)pn = Kp^, (4-5) 

where F and K are the polytropic exponent and polytropic constant, respectively. 
Up to the 2PN order, the four velocity is given by Eq.(3-35)E3)'cP'B). Since we need 
up to 3PN order to obtain the 2PN equations of motion, we derive it here. Using 
Eq.(3-35), we can calculate (au^)'^ up to 3PN order as 

+hijv'v^ + 2(^)(3'v' + (^\3)(3^v^ + 4(4)^ + - ^^)''^ + ^^^^ + ^'^l 
+0{c-'), (4-6) 
where we used 7*-^ = Sij — c~^hij + 0(c~^). Thus, we obtain up to the 3PN order 



as 



+^{-(6)« + I ((3)/5'(3)/?^' + h^.v'v^'^ + ^,)P^v^ + 5^s)(3^v^U - 2UX 

+0(c-^). (4-7) 

Substituting PN expansions of metric and matter variables into the Einstein 
equation, and using the polytropic EOS, we find that the metric variables obey the 



following Poisson equations I 



)• 



AU = -Attp, (4-8) 
AX = A-Kp (^v'^ + 2U + {Sr - 2)77 j , (4-9) 
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z\(4)V' = -27Tp[v'^ + n + 



Z^(5)/3^ = 167TP 



i I „,2 



V V 



+ 2t/ + ri7 +(3)/3' 



(4-10) 
(4-11) 

4f^,4(3)/5^- + (3)/3fi-|^^.-(3)/?,l. 



(4-12) 



Ah,, = {uU^.j - Ui.UAU - 3UiUj + 6ijU,kU,k^ - 16tt(^pvV - ^6ijpv^^ 



(3)/?:, + {3)/?,i - ^^^H3)P,k - 2 + 2(4)V'),i, - :,Si,A{X + 2(4) V) 



(4-13) 



zl(g)a = Anp 



2v^ + 2v^(5U + rn) + (3r - 2)i7;7 + 4(4) + X + 4(3)/3*v^ 



- h^JUij - -UUiUi + C/,z(2(4)V^ - X)^i 



+ 9(3)/5M(3)/5:.+(3)/?i-o'^^.-(3)/3,fc 



(4-14) 



Here, we consider the uniformly rotating fluid around z-axis with the angular 
velocity 12, i.e., 

v' = eijk^^x^ = {-yQ, xQ, 0), (4-15) 

where we choose fi^ = (0, 0, Q) and eijk is the completely anti-symmetric unit tensor. 
In this case, the following relations hold; 



d 



d 



dt dip 



d 



d 



dt dip 



d 



dt dip 



(4-16) 



where Q, Qi and Qij are arbitrary scalars, vectors, and tensors, respectively. Then, 
the conservation law (2-18) can be integrated asH) 



dP 



In n° + C7, 



pc2 + pn + P 

where C is a constant. For the polytropic EOS, Eq.( [4-17 ) becomes 



In 



1 + 



c^{r-iy 

Using Eq.(4-7), the 2PN approximation of Eq.( ^T^ ) is written as 



(4-17) 



(4-18) 



+ ^ ( -(6)« + ^(3)/5'(3)/?' + ^{3)P'v'U - ^ + {3)P'v'v^ + 2(4)V't^' 
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+C, (4-19) 



where H = rKp^-^/{r - 1), = R^f2^ and = + y^. Note that Eql4-19) 
can be also obtained from the 2PN Euler equation hke in the first PN caseE3). If 
we solve the coupled equations (4-8)-(4-14) and (4-19), we can obtain equilibrium 
configurations of the non-axisymmetric uniformly rotating body. 

In the above Poisson equations for metric variables, the source terms in the 
Poisson equations for (3)/3*, (5)/?*, and hij fall off slowly as r — > 00 because these 
terms behave as 0(r^^) at r — > 00. These Poisson equations do not take convenient 
forms when we try to solve them numerically as the boundary value problem. Hence 
in the following, we rewrite them into convenient forms. 

As for hij, first of all, we split the equation into three parts as& 

,{^) _ rrf r/- . . _ Ix. . ATT\ _ 'iTT .TT . X. TT . TT . = _A^q(^) 



Ah)^> = U[U,,, - -6,jAU ) - SUiUj + 6i,UkUk = -4^5>r' (4-20) 
Ah\p = -IGttLvV - y.jpv^ ] , (4-21) 



- 2(^(X + 2(4)^),,, - ^5,jA{X + 2(4) V^)) . (4-22) 

The equation for has a compact source, and also the source term of /i^j^^ behaves 

as 0{r~^) at r — > 00, so that Poisson equations for them are solved easily as the 

(G) 

boundary value problem. On the other hand, the source term of h\j behaves as 

0(r-3) at r ^ 00, so that it seems troublesome to solve the equation as the boundary 

(G) 

value problem. In order to solve the equation for h - as the boundary value problem, 
we had better rewrite the equation into useful forms. As shown by Asada, Shibata 
and FutamasecP, Eq.(4-22) is integrated to give 

/iSf =2^ J{pv^y\^-y\d'y + 2^ J (pv^yi^ - y\d'y + 5^, J - y\d'y 



12 dx'dxi 7 dx^dxi J V 2 



2 



pv 



2 + 3P - pU/2) 



I rr^^ -d'y- (4-23) 



Using the relations 



P = -ipv% + Oic-'), 

e = 0, 

v'x' = 0, (4-24) 
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Eq.(4-23) is rewritten as 



(G) 



7 



1 

+ 2 



|l(-'<5"'-<3f)+7^(--<3"'-'3 



d 



where 



-A-Kpv, 
AQiy = -Attpv'x^, 

AQ^^^ = -A'kI^pv'^ + 3P - \pU^, 



AQ 



(I) 



-47r( pv^ + 3P--pU ]x' 



-6^jQ^'\ (4-25) 



(4-26) 
(4-27) 

(4-28) 
(4-29) 



(G) 

Therefore, h\j can be deduced from variables which satisfy the Poisson equations 
with compact sources. 

The source terms in the Poisson equations (4- 11) and (4T2) for (3)/?* and (5)/3* 
also fall off slowly. However, we can rewrite them asa) 



(5)/?' = -4(5)^i - :T 2x^4)^ - 77i , 



where 



Aqt = —4:TTpx\ 



(4-30) 
(4-31) 

(4-32) 



^(5) Pi 



-Anp 
1 



i I „,2 



V V 



+ 2U + rn + (3)/3' 



+ f^4(3)AU(3)/?:-o'^u-(3)/?rfc 



{UU), - -((3)/3'^,^),^, 



Arji = -Anplv'^ + n + -U]x 



(4-33) 
(4-34) 



Thus (3)/3* and (5)/?* can be obtained by solving the Poisson equations (4-32)-(4-34) 
whose source terms fall off fast enough, 0(r~^), for numerical calculation. 

For later convenience, using the relation (3)-Pi = ^izkQk^ ^-nd Eq.( [4-16| ), (3)/5* and 
(5)/3* may be written in the form with explicit 17 dependence as 



(3)13' = ^ 



^(3)f3\ 



(4-35) 
(4-36) 
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where 

„fe I „,2 



j^x'iv' + 2U + rn\+^3)(3' 



+ ^>.({3)/3^ + (3)/3j-^5i,(3)At) 



+ l{UU^),-]i(3)^''Uk),. (4-37) 



4 

4.2. Basic equations appropriate for numerical approach 

Although equihbrium configurations can be formally obtained by solving Eq.(4-19) 
as well as metric potentials, U, X, (4)V', (6)Ck, (3)/5') (5)/3' and hij, they do not take 
convenient forms for numerical calculation. Thus, we here change Eq.(4-19) into 
forms appropriate to obtain numerically equilibrium configurations. 

In numerical calculation, the standard method to obtain equilibrium configura- 
tions is as follows E3^'E3); 

(1) We give a trial density configuration for p. 

(2) We solve the Poisson equations. 

(3) Using Eq.(4-19), we give a new density configuration. 

These procedures are repeated until a sufficient convergence is achieved. Here, at (3), 
we need to specify unknown constants, 17 and C. In standard numerical methods 
LJ^'IZHP, these are calculated during iteration by fixing densities at two points; i.e., 
if we put pi at xi and p2 at X2 into Eq.(4-19), we obtain two equations for 17 and 
C. Solving these two equations give us 17 and C. However, the procedure is not 
so simple in the PN case: 17 is included in the source of the Poisson equations for 

the variables such as X, (4)^"; (6)Q^' Vi^ (5) A) ^If^ ^IJ^' Q*'^^ Qi^'^ ■ Thus, if we 
use Eq.(4-19) as it is, equations for 17 and C become implicit equations for 17. In 
such a situation, the convergence to a solution is very slow. Therefore, we transform 
those equations into other forms in which the potentials as well as Eq.(4-19) become 
explicit polynomial equations in 17. 

First of all, we define q2, q2i, Qi, Qu, Qe and qij which satisfy 

Aq2 = -iirpR^, (4-38) 

Aq2^ = -iirpR^x', (4-39) 

Aq4 = -AnpR^, (4-40) 

Aqu = -47rpU, (4-41) 

Aqe = -Airpn, (4-42) 

Aqij = -Ajrpx'x^. (4-43) 

Then, X, (4)^/', Q^^\ Q'i \ Vi^ (5)A, Qij \ and are written as 

X = -2q2n^ -2qu-{3r-2)qe, (4-44) 
(4)i^ = ^(^q2^^ + qe + ^qu^, (4-45) 

= q^n^ + 3(r - l)qe - ^qu = q2^2^ + Q\^\ (4-46) 
Qf) = ^2,f22+Q(p, (4-47) 
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Vi = Q2i^ + mil 

Pi = ^izk<l2k^'^ + (b^Poii 



<3iP = ^iziqij^, 



where Qq^, ??oi and (5)^0? satisfy 

Z\Qf,? = -47r(3P - ^pU^' = -47rp(3(r - l)n - ^ 



U \x 



Am = -47rp{n + -u]x\ 



(4-48) 
(4-49) 

(4-50) 
(4-51) 



(4-52) 
(4-53) 



ei,kx''(2U + rnj+^2)P' 



+ ^{UU^),i-^{(3)P''Uk),i^-^7rS: 



(4-54) 



(G) 

Note that (5)/?* and /i^^ are the cubic and quadratic equations in i7, respectively, as 



where 



-4(5)Poi + UxHqe + lqu 



Voi,ifi 



{5)P'^^^ =-'^^izkq2k + l 



X q2,,f - q2i,ip 



(4-55) 
(4-56) 

(4-57) 
(4-58) 



[A) 



d_ 



"41'^ ^j'zfc^fc.V ^izkqk,<p ^izkqkj,ifi zkqki,ip ) ~l~ '^ij'^ ^kziqi 



+ 8^"i ^ I ^^^kziqi,<p - ekziqij,<fi ) + ^ I x'ekziqi,^ - ekziqii,^ 



d_ 

dx^ 



1 k\ d 
I ox 

Finally, we write ^Q•^a as 

(6)" = (6)"o + (6)«2^?^ - 2^4/2^, 
where (6)Q;o and (6)0:2 satisfy 



. (4-59) 



(4-60) 



Z\(6)ao = 47rp 
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= -4^5("«\ (4-61) 

+3q2,iUi + ^(3)/3j,- (^(3)/5^- + (3)/3fi - ^5ii(3)/3! 
= -4^5(^2)^ (4.g2) 

Using the above quantities, Eq.(4-19) is rewritten as 
where 

^ = ^ + ^ (29" + (3^ - 2)'?e) + ^{-(6)«o - ^ + f^(2'7« + (3r - 2)(7e) |, 
^ = ^ + ^ (2i?'C/ + 2(72 + (3)/3^) + ^{-(6)«2 + ^(3)/3\3)/3* + 4(3)/?^^ 

+(3r - + ^g^R' + ^[/^i?^ + 2g,C/ + + ^ + h^^J^^ 

D = ^ + ^{2(?4 + (3)/3^i?' + ^92i2' + 2?7i?4 + (5)^^) + ^ (/ij,^) + ^R\ri^Y 

(4-64) 

Note that in the above, we use the following relations which hold for arbitrary vector 
Qi and symmetric tensor Qjj, 

Qip — yQx ~l~ xQy, 
Qipip — y Qxx '^xyQxy ~i~ ^ Qyyi 
R^Qrr = x'^Qxx + 'i.xyQxy + y'^Qyy. (4-65) 

We also note that source terms of Poisson equations for variables which appear in 
A, B and D do not depend on Q explicitly. Thus, Eq.( [4-63| ) takes the desired form 
for numerical calculation. 

In this formalism, we need to solve 29 Poisson equations for U , Qx, Qy, Qz, {5)Pox, 

P n n d(^) n(^) n(^) n n n n n n h^'^^ h^'^^ h^'^^ h^'^^ h^'^^ 
{5)^0y^ VOx, VOy^ Vox ' '^Oy ' ^Oz ' ^2, q2x, q2y, q2z, Qu, Qe, iLxx , llxy , Hxz , Ityy , llyz , 

Qxx^ Qxy, Qxz, Qyz, (6)0^0) (6) and 54 . In Table 1, we show the list of the Poisson 
equations to be solved. In Table 2, we also summarize what variables are needed 
to calculate the metric variables U, X, (4)^^, {e)Oi, (3)/?% {5)/?*, hf^\ ^if^ and 

h[^\ Note that we do not need (^5)Poz, f/oz, and Qzz because they do not appear in 
any equation. Also, we do not have to solve the Poisson equations for /li^'' and qyy 
because they can be calculated from /li^^ = —h^xx — h^yy and qyy = q2 — Qxx- 
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In order to derive U, qi, q2, q2i, Qa, Qe and qij, we do not need any other potential 
because only matter variables appear in the source terms of their Poisson equations. 
On the other hand, for g^, QqP' Voi and hl^\ we need the Newtonian potential U, 
and for (5)-Poi) {6)'^o and (6)0:2, we need the Newtonian as well as PN potentials. 
Thus, U, qi, q2, q2i, Qi, Qe and qij must be solved first, and then QoP, Voi, 
(5)Poi and (6)02 should be solved. (6)Q;o must be solved after we obtain q^ because 
its Poisson equation involves qu in the source term. In Table 1 we also list potentials 
which are included in the source terms of the Poisson equations for other potentials. 

The configuration which we are most interested in and would like to obtain 
is the equilibrium state for BNS's of equal mass. Hence, we show the boundary 
condition at r — ^ 00 for this problem. When we consider equilibrium configurations 
for BNS's where the center of mass for each NS is on the x-axis, boundary conditions 
for potentials at r — > 00 become 



U = - I p(fx + 0{r-^), 
r J 

Q2 = lJ pR^d^x + 0(r-3), 
(le = \j pnd'^x + Oir-'^), 
qu = - [ pUd^x + 0{r-^), 



^- = % j px^d^x + 0{r-% 

^y = % I py^d^x + 0{r-^), 

^- = 5/ Pz^d^x + Oir-^), 
54 = ly" pR'^d^x + Oir-^), (4-66) 



(5)^0. = SPxd'x + ^ I SPyd^x + 0{r-^), 

{5)Poy = 5 / ^P^d'^ + S^yd'x + 0(r-3), (4-67) 

qS = ^Jp^' - 1)^ - \uy^x + 0(r-4), '?2. = 5 / pR^x^d^x + 0{r-% 

Qoy = 5 / Py'i^ir - 1)77 - ^Uy^x + 0{r-% <l2y = ^J pRVd'x + 0{r-% 

QoJ = ^ J pz'(3{r-l)n-^Uy^x + 0{r-^), q2z = ^ J pR'z^d^x + 0{r-^), 

(4-69) 

h(^J = 1 1 Sg)d'x + 0{r-% 4? = ^ / S^hyd^x + 0{r-% 
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h'^J = ^ I Sl,%zd'x + 0{r-% (4-70) 



1 f 3n^ny f 

1^^ = - j px^d^x + 0(r-3), q^y = —3— / px^y^(fx + ©(r"^), 

px^z^d?x + 0(r-5), gj,, = / py-^z^d^x + 0(r-5), (4-71) 



(6) 



ao = ^J S^'^°U^x + 0(r-3), (6)02 = ^J S^'^^^cfx + 0{r-^), (4-72) 



where n' = x'/r. We note that sP 0{r-^), S^f^ 0{r-^), ^("o) ^ 0{r-^) 
and S'^"^) — > 0(r~^) as r — > 00, so that all the above integrals are well defined. 

We would like to emphasize that from the 2PN order, the tensor part of the 
3- metric, 7ij, cannot Jie neglected even if we ignore gravitational waves. Recently, 



Wilson and Mathewsll^j), Wilson, Mathews and Marronettit21l' presented numerical 
equilibrium configurations of binary neutron stars using a semi-relativistic approx- 
imation, in which they assume the spatially conformal flat metric as the spatial 
3-metric, i.e., ^ij = 6ij. Thus, in their method, a 2PN term, hij, was completely 
neglected. However, it should be noted that this tensor potential plays an important 
role at the 2PN order: This is because they appear in the equations to determine 
equilibrium configurations as shown in previous sections and they also contribute to 
the total energy and angular momentum of systems. This means that if we performed 
the stability analysis ignoring the tensor potentials, we might reach an incorrect con- 
clusion. If we hope to obtain a general relativistic equilibrium configuration of binary 
neutron stars with a better accuracy (say less than 1%), we should take into account 
the tensor part of the 3-metric. 

In our formalism, we extract terms depending on the angular velocity 17 from 
the integrated Euler equation and Poisson equations for potentials, and rewrite the 
integrated Euler equation as an explicit equation in 17. This reduction will improve 
the convergence in numerical iteration procedure. As a result, the number of Poisson 
equations we need to solve in each step of iteration reaches 29. However, source terms 
of the Poisson equations decrease rapidly enough, at worst 0{r~^), in the region far 
from the source, so that we can solve accurately these equations as the boundary 
value problem like in the case of the first PN calculations. 

The formalism presented here will be useful to obtain equilibrium configurations 
for synchronized BNS's or the Jacobi ellipsoid. In this context, Bonazzola, Frieben 
and Gourgoulhon (1996) obtained an approximate nonaxisymmetric neutron star 
by perturbing a stationary axisymmetric configuration. However, they do not solve 
the exact 3D Einstein's equation. Thus, it is interesting to examine their conclu- 
sion on the transition between equilibrium configurations, which are approximately 
ellipsoids, by our PN methods. 
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§5. Gravitational Waves from post-Newtonian sources 

In the post-Minkowskian approximation, the background geometry is the Minkowski 
spacetime where linearized gravitational waves propagate EiP' The corrections 

to propagation of gravitational waves can be taken into account if one performs the 
post-Minkowskian approximation up to higher orders. In fact, Blanchet and Damour 
obtained the tail term of gravitational waves as the integral over the past history of 
the source O'EJ). They introduced an unphysical complex parameter B, and used 
the analytic continuation as a mathematical device in order to evaluate the so-called 
log term in the tail contribution. The method is powerful, but it is not easy to see 
the origin of the tail term. Will and Wiseman (1996) have also obtained the tail 
term from the mass qiiadrupole monient, by improving the Epstein- Wagoner for- 
malism eIP. NakamuraEy) and SchaferEy-* have obtained the quadrupole energy loss 
formula including the contribution from the tail term by studying the wave propa- 
gation in the Coulomb type potential. This indicates clearly that the origin of the 
tail term is due to the Coulomb type potential generated by the mass of the source. 
In the following treatment we calculate the waveform from the slow motion source 
and show explicitly how the tail term originates from the difference between the flat 
light cone and the true one which is due to the mass of the source GM/(? as the 
lowest order correction. 

5.1. Gravitational waves in Coulomb type potential 

We wish to clarify that the main part of the tail term is due to the propagation of 
gravitational waves on the light cone which deviates slightly from the flat light cone 
owing to the mass of the source. For this purpose, we shall work in the harmonic 
coordinate since the deviation may be easily seen in the reduced Einstein's equation 
in this coordinates. 

_ = -16^0/^^ + h^^'Ji"!^^^. (5-1) 

Since we look for 1/r part of the solution, the spatial derivative is not relevant in 
the differential operator so that it is more convenient for our purpose to transform 
Eq. (|5-l| ) into the following form 

(□ - h'^^dodo)h'"' = -167^5'^^ (5-2) 

where we defined S^'^ as 

= _ _ 2hVQ^ - h}^h^\ij). (5-3) 

Substituting the lowest order expression for derived in section 2, we finally obtain 
our basic equation to be solved. 



4M\ 92 



1 + + A 



r 



dt^ 



h"'' = (5-4) 



where the effective source f^'^ is defined as 

1 /^nn 4M\ 



= ^''-Y^ ( - — ) ^%o- (5-5) 
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At the lowest order, we obtain (4)T^!/ = {4)Ofj_iy. 

The solution for ( ^-41) may be written down by using the retarded Green's func- 
tion in the following form. 

h^"'{x)= J dx'G[^j\x,x')f^"'{x'). (5-6) 

The retarded Green's function is defined as satisfying the equation; 

aMG^^,\x,y) = 6\x-y), (5-7) 

with an appropriate boundary condition. We also defined the following symbol for 
the differential operator appearing in our basic equation. 

^2 



□ 



M 



4M\ 



r 



(5-8) 



where A is the Laplacian in the flat space. 

The Green's function Go for M = i.e. in the Minkowski spacetime is well 
known. We shall present the detailed derivation of the Green function in appendix 
B, and present here the result 

= J^e^'^' I dusgn{u)U+^^'^{x)^^''^'^*{x')e{r - r') 

Im '' ^ 

where * denotes the complex conjugate, sgn{uj) denotes a sign of uj and we deflned 

^+u;lmu) and llrSivlmU) as 



•^■^"""(x) = jHe— V'i^/(p;7)>lm, (5-10) 

y ZTT 

where p = ur and uf''^ is a spherical Coulomb function satisfying the outgoing wave 
condition at null infinity, and is also a spherical Coulomb function which is regular 
at the origin. 

Since we study an asymptotic form of the waves, we use Eqs.(S-7) and (-B-8) in 
appendix B to obtain the asymptotic form of the Green function as 

Im 

+0(r"2) for r^oo. (5-11) 

In addition, for slow motion sources, we evaluate the asymptotic form of the Green 
function up to 0{Muj) as 
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+0{MW)y-''^^'-^~''Hujr'yYUn)Yr^{n') + 0{ 

27r(2/ + 1)!! 



du! 



1 + 2Mu;|i ^- ^ - + In 2M^ + |s5n(i 



+i{ln\uj\ + C) I + 0(m2cj2 



s=l 



+ 0{7 



(5-12) 



where C is Euler's number. In deriving the above expression we have used the 
fohowing expansion for c; and ai in Muj. 



ci 



1 + ttMuj + 0{M^uj^) 



(2/ + 1)!! 

ai = 2Muj{c - ^ ^) + 0{M^uj^). 
Now we apply the formula©'© 



s=l 



(5-13) 



w / dve'^'^'lnv + i / — e 
Jo Ji V 

to Eq.(5-12), then we obtain 



-'^sgn{uj) -i{\D.\uj\+C), (5-14) 



r^2vr(2Z + l)!! 



xmi-2M(-El+ln2M)| + 2M(/; 



xe-^-(*-'^-*')(u;r')'y,™(/2)li:;,(/?') + ©(r'^). (5-15) 

3 worthwhile to point out that ln2M in Eq.(5-15) can be removed by using the 
idom to time translation. 

Here, we assume the no incoming radiation condition on the initial hypersurface 
so that we may take 

■ (5-16) 



lim e-'^^^^-'-'^-^'^nv ^ Q. 

V— >oo 



Thus we can make the following replacement 



J —e'"^"^-iujj dve'"^" Inv = dve'"^" Inv J —. (5-17) 

Inserting Eq.( 5-17] ) into Eq.(5-15), we finally obtain the desired expression for the 

rptarrlpd Orppn's fiinrtinn 



' dt 



(5-17) 



retarded Green's function 



i-iy 



r-2vr(2Z + l)!! 



dio 



l + 2Mfy--ln2M^ — 
V ^ J dt 



s=l 
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+2M( / dve"^''lnv)^ 



ln.j^+0(M 



- part of 
r 

G„(...') + 2M^i:/<'"{i<2s) + S;}«»"" 



V, X, X 



+ 0(r-2), 



(5-18) 



where we defined the spherical harmonic expansion coefficient of the flat Green's 
function as follows. 



(5-19) 



As a result, we obtain the waveform generated by the linear part of the effective 
source which corresponds to (CI) by BlanchetQ'. 



^ 1=2 



xMpqL-2{t -r -v) 



M,,,.,it - r) + 2M^ / dv[ln{^) + ^ 7} 

^-3pY"'aL-2|ea6(p'5g)feL-2(* " f) 



+o(m2) 



+ 0(r-'), 



r — v) 



(5-20) 



where -Pijpq is the transverse and traceless projection tensor and parentheses denote 
symmetrization. MpqL-2 and SpqL-2 are the mass and current multipole moments 
generated by the full nonlinear effective source f^u- They take the same forms as in 
Thorne (1980). 

5.2. Comparison with the previous work 

By using the post-Minkowskian approximation, Blanchet obtained the radiative 
mass moment and radiative current moment asEP 

'2 )i/r2\ 



(5-21) 



and 



y«(n) = ^^(m) + 2GM dvSl^^\u - ^^)| + + ©(G^M^^ 



(5-22) 
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where their moments Ml and Sl do not contain the nonhnear contribution outside 
the matter, P is a constant with temporal dimension, and ki and k'i are defined as 

1-2 

'^' = E7 + 777^^rT7^^' (5-23) 



^ 1 , 2/^ + 5^ + 4 
s ' l{l + l){l + 2) 



and 



1 ; - 1 



In black hole perturbation, the same tail corrections with multipole moments induced 
by linear perturbations have been obtainedEZP. Compared with Eqs.(5-21)-( ^-24] ), 
Eq.(5-20) shows that log term and l/-^ in the radiative moments (5-21) and (5-22) 
originate from propagation on the slightly curved light cone determined by Eq. (|5-8| ). 
That is to say, log term is produced by propagation of gravitational waves in the 
Coulomb type potential in the external region. It is worthwhile to point out the 
following fact: Only the log term has a hereditary property expressed as the integral 
over the past history of the source, since the constants ki and k'i represent merely 
instantaneous parts after performing the integral under the assumption that the 
source approaches static as the past infinity. 

However, Eq.(5-20) does not apparently agree with Eqs.(5-21) and (5-22), be- 
cause the definition of the moments {Ml, Sl} and {Ml, Sl} are different. We shall 
show the equivalence between our expression and that of Blanchet (1995), by calcu- 
lating the contribution from nonlinear terms like M x Ml or M x Sl- It is noteworthy 
that the luminosity of gravitational waves obtained by the present approach agrees 
at the tail term i.e. 0{c~^) with that by the post-Minkowskian approximation. 

^ (/ + !)(/ + 2) 1 
^ = h (^-1)^ /!(2^ + !)!!< > 
^4/(/ + 2) 1 
+ § (1-1) (i + 1)1(2/ + i);i < >■ f^'^^) 

5.3. Contributions from the nonlinear sources 

In order to evaluate the waveform produced by the nonlinear sources in f^u, it 
is enough to use the flat Green's function 

Im 

(5-26) 

Since the nonlinear sources have the form of either M x Ml or M x Sl, we have 
to treat the following type of retarded integral 



□ "1 



1,A fF{t-r) 



^^{-)^Q 



i-r'Y 



(g + i)! 



Post-Newtonian Approximation 



49 



q , . A 1 
nio + - — — rf'i<a,nQ_i> )— ^ F [t - r) 



2q+l 



(5-27) 



The evaluation can be made by using the formula which will be proved in appendix 
C. 



□ 



-1 



F{t - r) 



lim , 



+ n + iy.r{-k + / + 3 + 2n - A) 



n=0 



n!(2/ + 2n + 2) 



VJt-r) + 0(r-2). 

r 



(5- 



Although some terms of nonlinear sources may produce disastrous divergence in 
Eq.(5-28), we expect these divergent parts cancel out in total, which is explicitly 
demonstrated in the appendix. 

Using the above formula we find 



r \ r 



2(c? + l) 



+ 1 



riic 



2q + 



Y<^j<ai"'Q-l> 



- % -r) + 0(r-2) 
r 



(5-29) 



which is same with the expression (C2) obtained by BlanchetB. 

Applying Eqs.(5-29) to the nonlinear source f^i,, we can evaluate its contribution 
to the waveform. Together with Eq.(5-20), we obtain the total waveform as 



7! nL-2^MpqL-2{t -r) + 2M^ j dv\\n\ — \ + ni 



2M 



xMpgL-2it -r-v)^ - j^naL-2^eab(pSg)bL-2it - r 



+0{M^) 



+ 0{r- 



(5-30) 



In this section, we have derived the formula (5-30) for gravitational waves in- 
cluding tail by using the formula (5-20) and (5-29). We would like to emphasize two 
points on the derivation: First, in deriving Eq.(5-20), spherical Coulomb functions 
are used, since we use the wave operator ( [5-8D which take account of the Coulomb- 
type potential M/r. As a consequence, ln(t'/2M) appears naturally in Eq.(5-30). 
This is in contrast with Blanchet and Damour's method, where an arbitrary con- 
stant with temporal dimension P appears in the form of ln(w/P). Our derivation 
shows that the main part of the tail, which needs the past history of the source 
only through Inv, is produced by propagation in the Coulomb-type potential. Our 
method might make it easy to clarify conditions for wave formula including tail. This 
remains as a future work. 
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The second point relates with physical application: At the starting point of 
our derivation, the Fourier representation in frequency space has been used. Such 
a representation seems to simplify the calculation of gravitational waveforms from 
compact binaries in the quasi-circular orbit, since such a system can be described by 
a characteristic frequency. Applications to physical systems will be also done in the 
future. 

§6. Conclusion 

We have discussed various aspects on the post-Newtonian approximation mainly 
based on our own work. After presenting the basic structure of the PN approxima- 
tion in the framework of Newtonian limit along a regular asymptotic Newtonian 
sequence, we reformulated it in the appropriate form for numerical approach. For 
this purpose we have adopted the (3+1) formalism in general relativity. Although 
we restricted ourselves within the transverse gauge in this paper, we can use any 
gauge condition and investigate its property relatively easily in the (3+1) formalism, 
compared with in the standard PN approximation performed so farEiP' 
Using the developed formalism, we have written down the hydrodynamic equation up 
to 2.5PN order. For the sake of an actual numerical simulation, we consider carefully 
methods to solve the various metric quantities, especially, the 2PN tensor potential 
(4)hij. We found it possible to solve them by using the numerical methods familiar 
in Newtonian gravity. Thus, the formalism discussed in this paper will be useful in 
numerical applications. As an example of the application, we have presented the 
formalism for constructing the equilibrium configuration of nonaxisymmetric uni- 
formly rotating fluid. We have also discussed the propagation of gravitational waves 
explicitly taking into account of the deviation of the light cone from that of the 
flat spacetime and obtained essentially the same result with that of Blanchet and 
Damour for the waveform from post-Newtonian systems. 
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Appendix A 

Conserved quantities for 2PN approximation of uniformly rotating fluid ■ 



The conserved quantities in the 2PN approximation will be useful to investigate 
the stability property of equilibrium solutions obtained in numerical calculations. 
Since the definitions of the conserved quantities are given in section 3, we shall 
present only the results in this appendix. 
(1) Conserved massB'; 

= / pj^x, (A-1) 



where 



P*= P 



1 fl 



1 /3 



ens 



13 



45, 



(2)ADM mass 

Madm = j Aipd^X = J PADAld^X, 



(A-2)- 
(A-3) 



where 
PADM = P 



1 + ^ v' + n + -u 



1 



+^ ( f * + 9v'u + rnv' + 5un + ^u^ + ^^^^pw ] + o(c-«) 



• (A-4) 



4 2 

(3) Total energy, which is calculated from Madm — Af* in the third PN order H^; 

^ = y pEd^x, (A-5) 



where 
PE = P 

+ 



1 



1 / 5 



„^ + 77 _ + + + rv'n + 2un - -W + -o)/?* 



\2 

1 fll 



+ 



(A-6) 



It is noteworthy that terms including (5)/?* cancel out in total. 



54 



H. Asada and T. Futamase 



(4)Total linear and angular momenta: In the case K^^ = 0, these are calculated from 
(Wald 1984) 



Pi = — lim (p KijU^dS 

1 

~ 8^ 



(A-7) 



where Jj = {pc^ + pll + P)avPui. Up to the 2PN order, the second term in the last 
line of Eq.(^-7) becomes 



1 
1 

16^ 



hjk,i{'i)l3%d^x 



lim <(> hju ii3\(i'^ n^dS = 



(A. 



where we use hjk —>■ 0{r ^) and (a)/?-' — > 0{r ^) at r — > oo, and the gauge condition 
hjk,k = 0- Thus, in the 2PN approximation. Pi becomes 



Pi = J Pid^x, 



(A-9) 



where 



Pi = P 



+(3)/3' (v"^ + 6U + rn) + v' f 2(3)/?^?;* + 10(4) V + ^rnu 



+^c/2 + rnv^ + wuv^ + v^-x]}+ o(c-^) 



The total angular momentum J becomes 

J = j V^d^x, 

where = -ypx + xpy. 

Appendix B 

Construction of Green's function 



(A-10) 



(A-ll) 



In this appendix we shall present a detailed derivation of Green's function for 
□m- The following procedure is similar to that by Thorne (1980) for the Minkowski 
background spacetimeE^P. 
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The Green's function satisfying Eq.( ^-7D can be constructed by using the homo- 
geneous solutions for the equation 



□m1^ = 0. 

The homogeneous solution for Eq.( i3-lD takes a form of 
where we defined 

p = tor. 

Then the radial function fi{p) = pfi{p) satisfies 



dp"" 



+ 1 + 



4Mu; /(/ + 1) 



P 



fi{p) = 0, 



(B-1) 

(B-2) 
(B-3) 

(B-4) 



so that Eq.([g^ is a solution for Eq. d-BTl) . Thus we can obtain homogeneous solu- 
tions for Eq.( [B-lD by choosing fi{p) as one of spherical Coulomb functions; u[^\p; 7), 
Fi(p;j) and Gi{p\^) with 7 = —2Mlo. Here, we adopted the following definition of 
the spherical Coulomb functionEZP as 



Fi{p-^) = cie'Pp^+^F{l + l + i-i\2l + 2\-2ip), 

{p;-f) = ±2ie^''''cie^'Pp^^'^Wi{l + l±i-i\2l + 2\^2ip), 



,(±) 



Gi{p-^)=^-{u\-^^e^'^^+u\ 



(B-5) 



where q and ai are defined as 



\r{i + i + i-i)\ 



" ^ {21 + 1)! 

ai = aigr{l + 1 + i^). 



(B-6) 



Here, F and Wi are the confluent hyper geometric function and the Whittaker's 
function respectively. These spherical Coulomb functions have asymptotic behavior 
as 

Fi ~ sin — 7 In 2/9 — — /vr -|- o"/ 



Gi ~ cos ( /9 — 7 In 2/9 — -l-n + ai 



^ gxp 



ibi ( p — 7 In 2/> — —Itt 



for r — > 00, 



(B-7) 



and 



Fi ~ Qp'+i 

o, ' 



p ' for r — > 0. 



(2/ + l)q 



(B- 
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In terms of spherical Coulomb functions introduced above, we can construct the 
Green's function in the standard way. 

Gg(x,x') = ^6^^^' / du;sgn{Lu)(^'''^"'{x)^^''^'^*{x')e{r - r') 

Im ■' ^ 

where * denotes the complex conjugate, sgn{uj) denotes a sign of uj and we defined 
^^^^"^{x) and <Z^^'^'"*(a;) as 




The retarded Green's function is obtained for e = +. 




(B-IO) 



Appendix C 

Derivation of the formula (5.28) 



We evaluate the asymptotic form of the retarded integral as 



□" 



d x'Gq{x, x') 



n 



r'k 



where we defined Fu, as 



—i J dioLoe ^'^^hi(ijjr)F(^hL j r'^dr' ji{ur')- 

for large r, (C-1) 

F^ = ^j dte'^'F{t). (C-2) 



Since the Bessel function is defined as 



we obtain formally 



POO 

dyMby)e-''yy^^-' = Y^ 



(C-3) 



roo 

/ dye~ 
Jo 



nWiu + n + 1) 
g (-)"(6/2)'^+2"r(/x + u + 2n) 

n 



ay yfx+v+2n-l 



■^^ n!a'*+'^+2"r(z/ + n + 1) 



(C-4) 



Putting a = —i, h = 1, jj, = —k + 5/2 — A and u = I + 1/2 in Eq.(C-4), we obtain 
where we introduced small A G C in order to avoid the pole of r{—k + Z + 3 + 2n) 
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for k > / + 3. Thus we obtain 



r''^dr'ji(u)r') 



ijj 



k-3 



k+i+3 oo^ + ; + 3 + 2n - A) 



2 2^+ 



-E 



n=Q 



22"n!r(/ + n+ I 



(C-5) 



From Eqs.(C-l), and (C-5), we obtain 



22-n!r(/ + n + |) 



+0(r-2). 

This equation is vahd for A; > 3. For 3 < k < I + 2, we obtain 

1-1 



(C-6) 



(/ + A;-2)! r 

for 3 < A; < / + 2, 



(C-7) 



since 



(C-8) 

This formula for 3 < < / + 2 has been derived by Blanchet and Damour (1988). 



Appendix D 

Derivation of (5.29) 



From Eq.(5-27), we shall derive Eq.(5-29) in this appendix. First, for the first 
term in the right hand side of Eq.(5-27), we use Eq.(C-7) to obtain 



^^«^lt-r) + 0(r-) 



1 

Uio^T^ F[t-r) 



2(g + l) r 

Next, the second term in the right hand side of Eq.(5-27) is calculated as 



(D-l) 



9-2 



(g + i)! 



-□ 



-1 



2q + l 



1 '^»<a,nQ_i> (g) 

F(t-r) + C'(r ), 



(D.2) 



and 



(-) 



_ (g + i)! 

2g + l .^i2Jj!((?-j)! 



' . 1 ('J-i) 
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{2q + l)ql 



X lim 



+ 



{2q)l 



r(-A) + r(-i-A) 



+ 0(r-2) 



Here we used Eq.(5-28), 



E 



n!(2g + 2n)! 



r{2n)+r{2n-l] 



2q(2q)\ 



and 



lim (r(-A) + r(-l - A)) = lim(-A)r(-l - A) = -1, 



which, in fact, means the cancellation of the poles. 
From Eqs.(5-27) and (D-l)-(D-3), we obtain 



□ 



-1 



2((? + l) 



F{t - r] 



+ 1 



n, 



i> 



1 H ^ 

- Fit -r) + 0{r' 



(D-S) 



(D.4) 



(D-5) 



(D.6) 



which equals to (C2) obtained by BlanchetLF. As for Eq.(D-6), we used unphysical 
but small parameter A in order to avoid the pole of the integrals. However, in the 
limit A ^ 0, we showed explicitly the cancellation of poles and obtained the finite 
values in total. 
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Table 1 

List of potentials to be solved (column 1), Poisson equations for them (column 2), 
and other potential variables which appear in the source term of the Poisson equation 
(column 3). Note that i and j run x,y,z. Also, note that we do not have to solve 

Voz, {5)Poz, Qyy, Qzz and h^^\ 



Pot. 


Eq. 


Needed pots. 




Pot. 


Eq. 


Needed pots. 


U 


(2.11) 


None 






(4.6) 


None 


Qi 


(3.14) 


None 






(4.15) 


U 




(4.1) 


None 




mi 


(4.16) 


U 


q2i 


(4.2) 


None 




(5)^0j 


(4.17) 


U, qi 


Q4 


(4.3) 


None 




(6)«0 


(4.21) 


'^■i Qey qui ii'ij 5 *Voi 


Qu 


(4.4) 


U 




(6) "2 


(4.22) 


U, 92, qi, 92 j, qij 


Qe 


(4.5) 


None 




"'^3 


(3.1) 


U 



Table 2 

Variables to be solved in order to obtain the original metric variables. 



Metric 


Variables to be solved 


see Eq. 


U 


U 


(2.11) 




qi, u 


(3.17) 


X 


92, qu, 9e 


(4.7) 




92, qu, Qe 


(4.8) 




(5)-Po«, VOi, Qu, qe 


(4.18) 




92i, 92 


(4.18) 


(6)" 


(6)«0, (6) "2, 94 


(4.20) 






(3.1) 


'hj 


qij,q2 


(4.14) 


'hj 


Qoi^, Qu, Qe 


(4.19) 


% 


qij, 92, 92i, 9i 


(4.19) 



